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CHAPTER I 



INTRODUCTION TO DISCRETE ORDINATES METHODS 
IN TRANSPORT THEORY 



Transport Theory 

The behavior of a nuclear reactor and the effective- 
ness of its shield are governed by the distribution in 
space, velocity, and time of the neutrons and photons in the 
system. The transport equation is a statement of conserva- 
tion for these neutral particles, the solution of which 
determines their distribution in phase space. Derivation 

of the transport equation is presented in most transport 
1-3 

theory texts . 

4 5 

A few special systems have exact solutions, ' but 
for more practical systems, numerical solutions of an 
approximate transport equation are sought. A few of the 

more useful transport approximations are spherical harmonics^ 

7 8 9 

Monte Carlo, discrete ordinates, and moments methods. 

The spherical harmonics method treats the anisotropy 
of the particle distribution and cross sections to various 
degrees of approximation. Its application to plane, spher- 
ical, and cylindrical geometries has been useful but for 
more complex geometries , the spherical harmonics method is 
so complicated that other methods are used. 
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The moments method, used in shielding studies, 
calculates particle transport rigorously but the method is 
limited to infinite homogeneous media. 

The discrete ordinates and Monte Carlo methods are 
the most nearly rigorous for multidimensional problems but 
the former is hampered by ray-effects^® and is currently 
limited to two-dimensional problems. The Monte Carlo method 
is the only useful technique for rigorous solution of 
three-dimensional problems. 

Development of Discrete Ordinate Method 

The use of discrete ordinates was first suggested by 

Wick in 1943^^ and developed for radiation transport in 

12 1 3 

stellar atmospheres by Chandrasekhar. ' Carlson intro- 
duced the angular "segmented" approximation for radiation 
transport calculations at Los Alamos in late 1952 and early 
1953. A complete history of development is given by 
Carlson^^ and a current bibliography is available.^® 
Chernick^^ gives a detailed historical review of reactor 
physics calculation methods from the Manhattan Project of 
the 1940 's to the present and describes the place of 
computer codes in reactor criticality and shielding appli- 
cations . 

Early methods have been refined to the present 
"diamond difference" scheme, described by Carlson,® which 
is incorporated in most codes. New approximations, 
such as MS^j , have been derived for special problems. 
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In 1967 Shreiner attempted to derive an Sfj differ- 
ence scheme, called VSj^, from a variational principle but 
produced no improvement over existing methods. More 
recently, effort has been expended to derive discrete 
ordinate approximations from space-angle synthesis tech- 
niques.^^ Kaplan^^ derives an S £ approximation in K,V 
geometry in this manner by using step functions in direction 
as trial functions to eliminate ray effects. Natelson'^'^ 
derives discrete ordinate approximations from a first-order 
variational principle in X,V geometry. For some unknown 
reason, this approach has not led to improved difference 
schemes over the standard diamond difference method. 

Applications of Sjij Method 

Transport theory, as opposed to one of the simpler 
approximation techniques such as diffusion theory, must be 
used when either very precise solutions are required or 
flux variations are large. As a rule of thumb, transport 
theory is used when the logarithmic gradient of the particle 
flow density is of the order of or larger than an inverse 
mean free path over significant portions of the system. 
Early application of by Lee^^ to the one-speed problem 
revealed that the Sq approximation gave critical radii 
within 0.3% of the values obtained by the exact solution of 
the neutron transport equation. Problems exist, however, 
for which the one-speed approximation is not sufficient. 
Mills applying to small fast critical assemblies such 
as Godiva, Jezebel, and Topsy, revealed that 24-group, Sg 
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calculations with a ?£ approximation to the scattering 

anisotropy should be sufficient to provide an accurate 

treatment of neutron transport in fast metal assemblies. 

More recently, the demand for tighter safety 

standards has awakened renewed interest in application of 

2 5 

transport theory Sf^ approximations. Protsik compared 
results of multigroup diffusion theory and S|^ calculations 
and found that significant error in loss-of-coolant reac- 
tivity calculations occurs using diffusion theory. He 
found that diffusion theory over-predicts the neutron leak- 
age and under-predicts the multiplication. 

The $ 1^1 approximation has also been used as a standard, 
for problems without an exact solution, by which the merits 
of other techniques are judged. solutions of the 

adjoint equation have been used in importance sampling tech- 
niques for Monte Carlo methods.^® 

Multigroup Equations 

The steady-state multigroup transport equations can 
be written^® 

Li|j(r,fi) - Sij;(r,n) + q(r,^) (1) 

where the vectors and q contain the unknown distributions 
and sources in each energy group, respectively. The opera- 
tors L and S take the following forms 

(Li{/(r,n))g = n‘V\pg(r,Q) + o (r)ipg (r,^) (2) 

(Si|»(r,n)) = J /dn' ^(r;n\g' ^ ,g) , (r ,5 ' ) (3) 

gl Q 

where the subscript g denotes the g^^ energy group. 
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The approximation discretizes the above equations 

A. 

at N discrete angular ordinates, Q yi, to produce the set 

L\p(r,n^) = (4) 

in which 

(Lij; (r,5 ^) ) g = ^ ^ 

and 

W ^ ^ 

g ' m=i 

where are weights associated with the angular quadrature 

set and the sum over m approximates the integration 

^ I 

over Q in (3). Equations (4), (5), and (6) describe the 

Sj^ approximation to the transport equation. 

Finite-difference techniques using the diamond dif- 

O 

ference approximation have been applied which approximate 
equation (4) to second order in all independent variables 
in Cartesian geometries and to order Ar /r in one- 
dimensional spherical geometry. Applying any method of 
spatial discretization, the Sf^ equations can be written in 
the matrix form 

L}l!_ = SjlJ_ + ^ (7) 

where the vectors \j^ and ^ now contain the unknown flux and 
inhomogeneous source, respectively, in each group and at all 
mesh points, space and angle. The matrices L and S are 
approximations to the operators L and S respectively. 
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Outer Iteration 

The scattering matrix S is split as follows 

S = ‘ (8) 

where S^,S^, and represent down scattering, self or 
in-group scattering, and up-scattering respectively. The 
iterative technique 

(9) 

represents the "outer" iteration where j is the iteration 
index. This matrix equation is never explicitly solved in 
transport codes. The matrix is dense and of such 

large dimension that it is not easily invertible on present 
day computers. Consequently another splitting is made so that 
an inner iteration can be conducted within each energy group. 



Inner Iteration 

The self or in-group scattering matrix is separated 
from the left hand side of (9) to operate on the inner 
iterate group flux. 






y + / , fe y + / 

^g'*'g " ^g 



(10) 



where the source to the group is 



In this inner iteration scheme, the first term on the right- 

/ + / 

hand-side of (11) can be evaluated even though \j^-^ 



is not 
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known for groups g and below, since elements of corre- 
sponding to these groups are zero. The matrices Lg and 

in (10) are restrictions of L and to the group. 

The vector 'Pg^^ contains components of corresponding to 

the group. The matrix Lg is triangular and thus easily 

invertible for typical difference schemes, so 
be obtained rapidly from ^^g* ^ ^ 

A converged inner iterate solution in each energy 
group completes one outer iteration. Consequently the 
outer iterate equation (9) is never solved but its solution 
is approximated by completing the inner iteration for each 
group a number of times, each with a source, £, calculated 
from previous outer iteration fluxes. For media which do 
not reproduce neutrons, only one outer iteration is required 
but for fissioning media, many outer iterations are required. 

The problems treated here employ the inner iteration 
of equation (10) for plane geometry. The specific equations 
will be developed later. 



Convergence of Sf^ Approximation 

Early versions of the discrete ordinates equations 
were used and results compared to those of other methods 
without regard to their formal convergence properties . 
Numerical experiments with problems representing physical 
systems showed that the early (circa late 1950 's) time- 
dependent schemes were stable. Not until 1960 did proof 
of convergence of the steady-state approximation appear 
when Keller^^ proved mean-square convergence of the Sfj 
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approximation to the analytic solution in plane geometry, 
one-speed, with isotropic scattering for values of 
0 ■' C < 1 where C is the mean number of secondary neutrons 
per collision. At the same time, Wendroff^^ proved point- 
wise convergence for the same problem for C < 1 and a 
weighted mean-square convergence for C = 1. Keller 
immediately thereafter proved pointwise convergence for the 
problem for values of C greater than unity. 

Not until 1968 did convergence proofs appear for 
multidimensional geometries. Madsen^^"^® proved pointwise 
convergence of the approximation to the time- independent 
one-speed angular segmented transport solution in x,y 
geometry with vacuum and periodic boundary conditions. In 
addition he proved that the spatially differenced approxi- 
mations using the central difference, first-order difference, 
and diamond difference schemes converge pointwise to the 

3 9 

angular S|y approximation solution in x,y geometry. Madsen 
has shown that the Sf^ approximation in x,y,z geometry con- 
verges pointwise to the exact solution under certain condi- 
tions . 

Acceleration of Convergence 

Current work is being directed toward improved tech- 
niques to accelerate convergence of the Sf^ algorithm. 
Clifford^O gives a lucid account of S|y running times for 
shielding problems using various generation digital com- 
puters. With present generation IBM/360 series machines, 
very large two-dimensional problems require 10 to 20 hours 
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running time. For deep penetration problems in media with 
a®/a^ near unity, Reed states that non-accelerated 
methods converge so slowly that they are practically use- 
less . 

Several methods have been developed to accelerate 

inner iteration convergence of the algorithm. Carlson 

and Bell^^ proposed a scale factor technique which is used 

in present standard codes like ANISN. It uses a system- 

wide neutron conservation principle and iterates until a 

o 42 

"false" source is arbitrarily small. Engle found that 
the scale factor technique accelerates convergence rapidly 
for absorbing media and near source regions for scattering 
media but is slow to converge at points many mean free 
paths away from the source in predominantly scattering media. 
For the latter problems, the scale factor rapidly approaches 
unity (even when the group flux solution is far from con- 
verged) and thereafter does little to accelerate conver- 
gence. Consequently, he proposed a separate scale factor 
for each space interval which exhibited some success in a 
one-dimensional 26-group problem. Clancy^^ proposed an 
outer iteration scaling technique which, when used with 
inner iteration scaling, accelerates convergence for prob- 
lems characterized by the presence of significant upscat- 
tering. 

The acceleration technique often used in one-dimen- 
sional problems is that of Chebychev, adapted to transport 
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approximations by Hageman.*^^ It is not as effective in 
two-dimensional geometries. 

The Synthetic Method, first proposed by Kopp^^ and 
applied by Gelbard^^ has been quite useful in accelerating 
convergence of the S|^ algorithm in one and two-dimensional 
problems. This method uses the solution of a diffusion 
theory approximation to accelerate convergence of the 
iterative S^j algorithm. 

Another effective method of accelerating convergence 
is called coarse mesh rebalancing. An early version was 
suggested by Wachspress^^ for diffusion theory codes. The 
method has been applied to a transport code, TWOTRAN.^® A 
variational rebalancing scheme was devised by Nakamura'^^ 
and applied to the one-dimensional Sfj algorithm. 

The most effective acceleration methods proposed to 
date are the synthetic and coarse mesh rebalance methods. 
When successful, they lead to tremendous reductions in com- 
puting time. Unfortunately, Reed discloses, there are 
problems for which the use of these acceleration techniques 
lead to an unstable algorithm. This failure occurs on 
those problems where it is most necessary to accelerate 
convergence, namely problems with an optically thick region 
with scattering ratio near unity. He displays model prob- 
lems for which the unaccelerated algorithm requires more 
than 2600 iterations to converge and the above two methods 
fail to converge. He also generates parameters which force 
both methods to converge. 
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A typical problem for which acceleration techniques 
are required is neutron propagation through large thick- 
nesses of iron which has a huge scattering resonance at 
about 25 kev. Devillers^® compares Sj^ codes running times 
(NIOBE, 20 minutes; ANISN, 8 minutes) to a Monte Carlo code 
(POKER, 13 minutes) for a one-dimensional multigroup prob- 
lem. Deep penetration problems in iron have also been 
solved by Sfj techniques to study the "window" effect for 
such materials. 

Current research is directed toward developing more 
accurate transport results. Increased emphasis on reactor 
safety is one factor responsible for this trend. But it is 
not sufficient to strive for accuracy alone; efficiency or 
cost must also be considered. Consequently, acceleration 
techniques play a prominent role in present computer codes . 
Any technique which contributes to improved efficiency will 
be incorporated to reduce the expensive computer computa- 
tional costs. 

Proposed Investigation 

This investigation formulates the angular flux form 
of the Sfj equations in plane geometry within a mathematical 
frame which allows an insight into the mechanics of the 
convergence process. This mathematical formulation leads 
to sufficiency conditions under which iterative convergence 
to an exact discretized solution is assured. These condi- 
tions include a domain in which some elements of the 
iteration matrix are allowed to be negative. 
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The infinity-norm of the iteration matrix is calcu- 
lated for a problem with a homogeneous isotropically scatter- 
ing medium. This allows the calculation of an improved con- 

CO 

vergence criterion. Reed states that the convergence 
criterion 



max 

■i 



1 - 



TTfe) 



< 0.0001 



where is the component of the scalar flux 

iterate is not, strictly speaking, a measure of the error in 
the final iterate , since slowly converging methods may meet 
this test without being close to the exact solution. It 
is, however, the convergence criterion actually used in 
most transport codes. This investigation proposes an 
improved convergence criterion which, along with the 
demonstrated convergence properties of the Sj^ algorithm, 
guarantees that the fractional iterative error, defined in 
Chapter II, is arbitrarily small. This is exhibited in 
Chapter II. 

A spatial transform method which accelerates inner 
iteration convergence for certain problems is developed 
in Chapter III. The method is applicable to problems for 
optically thick media in which the scattering ratio is near 
unity. The transform renders the angular discretized Sf^ 
equations invariant and places an acceleration parameter in 
each non- zero matrix element of the iteration matrix such 
that its norm is reduced. This results in accelerated con- 
vergence of the Sjy algorithm. 
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The reciprocal of the spatial transform is applied 
to an problem to reduce the discretization error. This 
variation of the transform method is applied to a predom- 
inantly absorbing medium problem in Chapter IV. 

An upper bound is found for the round-off error in 
the Sj^ iterative algorithm describing an inner iteration of 
a multigroup problem. The expression for this bound pro- 
vides insight into the reasons why iterative techniques 
often exhibit smaller roundoff error than competitive matrix 
inversion methods and how these errors propagate through the 
iterative process. This exposition is presented in Chapter 
V. 



CHAPTER II 



DISCRETE ORDINATES AS A METHOD OF 
SUCCESSIVE APPROXIMATIONS 

To analyze the convergence properties of the Sf^j 
algorithm, it is necessary to provide some theoretical back- 
ground and to show that the algorithm possesses, under 
certain restrictions, sufficient conditions for convergence. 

Theoretical Background 

The problem to be solved is 

( 1 ) 

where L is a finite dimensional square matrix, is the 
solution vector, and £_ is the inhomogeneous vector. All 
vector and matrix elements are real. For the applications 
considered here, L is non-singular so equation (1) has a 
unique solution, , for each given 

Now assume that L can be split as follows, 

L = D-E-S . (2) 

The choice of the splitting is arbitrary but the implication 
is that (D-E)"^ is an approximation to L~ ^ but (d-e) is much 
simpler to invert than L. For instance (D-E) may be chosen 
such that it is sparse compared to the dense L. 
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After splitting L, equation (1) becomes 
(D-E) tl-(D-E) ■ ^S]^* = £_ 



or 



^* = [I-(D-E)"^S]‘Md-E)"^£ 



where 



L'^ = [I-(D-E)’^S]"^ (D-E)"L 



Use will now be made of a convergence theorem from 
53 

Isaacson, which is repeated here. The notation has been 
modified for convenience to conform to that of the Sfj 
algorithm used later. 



Theorem 1 

The geometric series 

only if 



niD-E)"^S]"’ 



m=o 



converges if and 



I I (D-E)-'S| I < I . (3) 

If (D-E)'^S is convergent, then [T-(D-E)~^S] is non-singular 
and 

oo 

[I- (D-E)- ^S]-^ = . 

m=o 

If the splitting of L is such that expression (3) is 
satisfied, the theorem can be invoked to produce 

CO 

L"^ = ^ t (D-E)- ^S]"’(D-E)- ^ 
m=o 
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which is the Neumann series expansion of L'^ and 

00 

= n<D-E)"^5]"'(D-E)"^2. • 

m^o 

A consequence of the series expression for L~^ is that 

/- 1 1 Id-e) ■ 

Observe that | | L” ^ | | is bounded provided that | | (D-E)"^S| | 

< 1 and that the bound becomes tighter as | | (D-E)"^S| | 
approaches zero. 

Iterative Process 

On the basis of Theorem 1, Rall^^ constructs an 
iterative process, called the method of successive approxi- 
mations, which has 

= (D-E)‘^S^f^J + (D-E)'^2^ (6) 

where -i is an iteration index. This process gives a 
sequence of successive approximations which converge 

to the exact solution starting from any initial guess 
Equation (6) implies that 

^(^+/) ^ J^[(D-E)'^S]'”(D-E)‘^£ + [(D-E)'^S]'^‘^^^^^L (7) 

m= 0 

Observe that t^m since the first term on the 

-C 00 

right hand side of (7) approaches and the second term 
vanishes due to condition (3) . 
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Error Estimation 

Using equations (4) and (7) , it follows that the 
absolute error, defined by 

e (8) 



is bounded by 

||e<-^l|| < J ll P - - .£- ) .l ' .Si r^- . |||D-E|-',|| ,9, 

;-| I (D-E)"'S| I 

for ^ = 0 . 

This error expression is of little practical value except to 
illustrate the idea of efficiency. It is desirable to 
obtain approximations of a given accuracy with the fewest 
number of iterations. Obviously, more iterations are 
required to achieve a given accuracy as | | (D-E ) ^S| | 

approaches unity. Indeed the error is unbounded for 
||(D-E)"^S|| = /. But this condition violates expression 
(3) which guarantees that the series converges to L~K If 
I I (D-e)”^S| \ ^ 1 , the convergence of the iterative process 
is not guaranteed. 

Formulation of Sfj Algorithm 

The equations will be derived for an energy group 
in a multigroup problem applying the following assumptions 
to the transport equation: 

1) steady state 

2) infinite plane geometry with one spatial dimension 
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3 ) homogeneous medium 

4) isotropic scattering in the Laboratory Coordinate 
System 

5) no reentrant current at the slab boundaries. 

The problem may be described by 



/ 



^ 9X 


/ 0^ r 

+ 0 = 2 i 


i//(x,y' )dy' + q |x,y) 


(10) 




- / 






with boundary 


conditions 








'l^(o,]i) = 0 


for y > 0 


(11a) 




(i, y ) = 0 


for y < 0 


(11b) 



where X = 0 at the left face of the slab and X = L at the 

right face of the slab. Equation (10) is the balance 

equation for the angular flux in a particular energy group. 

The source term q(X,y) represents neutrons scattered from 

other energy groups and born within the energy group, either 

from fission or external sources. 

The angular and spatial discretization of equation 

(10) follows the procedure described by Bell.^^ Alternate 

5 6 

formulations by Mynatt, who applies numerical difference 

O 

methods to the analytic balance equation, and Carlson, who 
immediately writes a discretized form of the balance 
equation with unknown coefficients and evaluates them by 
applying neutron conservation principles, produce the same 
difference equations in plane geometry. Equation (10) is 
discretized in direction cosines by approximating it at each 
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of the W values of yy in the quadrature set selected. This 
gives the N coupled set of equations 

d\l> j{x) 4 W 

y. — ^ — + o^ipflx] = ^ + q - lx) 

J n=? J 

for j = 1,2,***, W where ' '^J 

for y = 1,2,**‘, W/2. 

A spatial mesh is selected with R equi-spaced incre- 
ments providing R+ / space points, including the boundaries, 

at which ij^ (x) is to be approximated. Each increment. A, is 
i 

defined by 

A = X, - X, . (14a) 

k+ 1 k 



The x,y mesh is described in Figure 1. 

The quantity ipj(x) in equation (12) is approximated 
at each spatial midpoint of the mesh 





X, , 

^fe+//2 " 2 


(14b) 


by 












(14c) 


and the 


derivative term is approximated by the 


central dif- 


ference 


dip: ( X ) 
cCx~ 


*k*Ki ■ *k.j 
“ ^ 


(15) 




^k+112 
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FIGURE 2.1 

SPACE-DIRECTION COSINE MESH 




Direction set properties; 

1) • -y y • ; y . > 0. 

■ N 

I W ■ /, U)jy ■ 0) . , y • J,2, • ♦ ♦ ,y ; Wf > 0. 

nmj TT+J ^ 



2 ) 
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This approximation of ^jix.) j ^ type 

5 1 

of error called discretization error. Isaacson shows 

2 7 

that this approximation is accurate to order A , O(A^). A 



desirable strategy is to choose A small enough that dis- 
cretization error is reasonably small. This choice is based 
on a guess of how the solution is expected to vary across 
the slab. 

The above substitutions applied to equation (12) 
produce the following set of equations . 

For ]} j > 0 




For yy < (? 









k = R,R- In these equations 



2|yy| + Ao"^ 



(16c) 



72 : 



and 






TE 



(16d) 



and 







(16e) 



Equations (16) could be arranged in the matrix form 
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and the solution found by inverting L. This procedure 
is not useful, however, because L is usually dense and is 
of such large dimension for practical reactor problems that 
present day computers do not have sufficient storage 
capacity to invert such systems. Instead of the above pro- 
cedure, an iterative technique is used for reactor problems 
which solves (16) in the following manner. 



For y • > 0, 
J 



J k+1,J 



j^k,j 



/ w 



(^) 



— yw (ii > 0 , . , 
4 ^k,n' 



(17a) 



and for ]Sj < 0, 



j I (•< + / ) I (-C+ / ) 

d.ii; - e-ip, . . 

J k,j 






where -c is an iteration index. These equations are the 
basic equations describing the algorithm. 

The iterative solution of this system, , approx- 

imates the exact matrix solution ^*. The error in this 
approximation is called iteration error. A practical 
strategy is to iterate until this error is arbitrarily small. 



Calculation Sequence 

The convention of performing the calculation 
spatially in the direction of neutron flow, as pointed out 

C C 

by Bell, has been observed in equations (16) and (17) . 
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Consequently , iterative and round-off errors are attenuated 
rather than amplified as the calculation proceeds. This is 
observed from the S^j algorithm in its machine form. That 
is, for Uj > 0 



(/C+ / 

fe+;,y 




(i+n 




w 



n=1 fe+J/2,n 



■k+1 12 ,j) 



(18a) 



and for u • < 



0 



k, 




(i+ 1 ) 
fe+ I , j 




hfk^’/2 



U) 






(18b) 



e/ 

The quantities are less than unity for all j so any 
numerical errors in a particular angular flux value will be 
attenuated in calculating the neighboring angular flux value. 
Amplification of these errors would have occurred had the 
convention not been observed. 



Iterative Matrix Formulation of the S^j Algorithm 
Equations (18) with boundary conditions 





for yy 


> 0 


(18c) 




for yy 


< 0 


(18d) 



and an initial guess = 0 for all k,j is the machine 

algorithm for the problem described by equations (10) and 
( 11 ) . 

Equations (17) are precisely the same S^j algorithm 
expressed in slightly different form. If these equations 
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for an R+/ by W mesh (including boundary points) arc written 
so that the solution set is arranged in the following vector 
order 



'P2,2 

'^R+I ,2 

ih * 

^R+I,M/2 
’^'R,N/2+I 
’^R- I,W/2+/ 

,W/2+/ 
'^'R,N/2 + 2 
’l'R-7,W/2 + 2 

,H/2 + 2 




(19) 



which has dimension R*N, then equations (17) describe the 
matrix system 



+ 2. . 






( 20 ) 
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The matrix (D-E) is block diagonal of the form 




( 21 ) 



where each block is an R by R submatrix. R is the number 
of spatial increments in the mesh. There are N such blocks. 
Each block By has the lower bi-diagonal form 




-®y 







( 22 ) 



where / is the index belonging to py in the direction set 
from which Sy and dy are calculated by equations (16c) and 
(16d). 
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The S matrix also has a regular block structure of 
the form 





^BP^ 


• • • 




r’ 


• ♦ ♦ 


^«W-I 






^BP2 








^CPw 


• 




♦ • 

• • 

• « 




• ♦ 

• • 




^BP2 


• ♦ ♦ 

• « • 


^^^N/2 




• • '• 

♦ ♦ • 






^CPj 


^CP^ 


^^^N/2 






"8^ 


* * 

♦ • 

• • 




• • 

• « 

• • 




♦ « 

• • 

• « 


^CP, 


^CP^ 


• • ♦ 

• • • 


^CP^/2 




♦ • • 




^BVn 


c 

^CPj 


^CP^ 


^CV 

^^NI2 


J+l 


^sv, 





where the notation and Sg^ distinguishes the form of the 

block submatrix and the subscript denotes the index belonging 
to all elements in a particular block. The S matrix has 
certain block regularities listed below: 

1) The same index appears in each column of blocks. 

2) Each block is an R by R submatrix of S. 

3) There are W blocks in each row and column of 

blocks . 

4) The upper right and lower left quadrants of S are 
composed of blocks, each with form Sq-q. 

5) The upper left and lower right quadrants of S are 
composed of blocks, each with the form Sjgp. 
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The block Sgp_ form is 

/C 



6 . 4 . 






O 



o 



^ . 4 . 

-c -c 



( 24 ) 



where each matrix element has the same value 






( 25 ) 



and is the weight in the quadrature set. Observe 

that 5gp. is lower bi-diagonal. 

The block S^p. form is 



O 






'A, 









O 



( 26 ) 
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which will be called upper cross bidiagonal. The matrix 
elements are given by (25) . 

Inverting the Matrix (D-E) 

Due to the simple block form of (D-E) and the fact 
that each block is lower bidiagonal, see (21) and 22) , its 
inverse is of a simple form. In fact (D-E ) ^ is of the 

block form 




(27) 



where each block. By, of (21) is inverted individually and 
placed in the same position in (D-E)”^ that it held in 
(D-E) . 

Each block can easily be calculated from (22) 

57 - J 

using the method described in Wylie. By has the form 
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l^th 



/LOW 



J 









2 




1 

d- 


3 


4 










efe- 3 


d^ 

■c 




J-2 


eK- 1 






-A 




-c 


P 


~T-J 


T-2 






dX 



o 



i 



= 8 



-I 



(28) 



— 

.2 



which is lower triangular with elements given by (16c) and 
(16d) . 

Since (D-E)"^ and S are now known, equation (20) can 
be operated on by (D-E)”^S to obtain 



^(-<-*1) = (D-E)'^S + (D-E) ^ 2^ . (29) 

The matrix (D-E)"^S is called the iteration matrix 
and can be determined from (23) , (24) , (26) , (27) , and (28) . 

(D-E) can be formed by multiplying (27) onto (23) block 
by block: That is, due to the block diagonal form of 

(D-E)”^, (D-E)”^S is formed by matrix multiplying 8^^ onto 

each block in the row of blocks in S. That is. 
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^ J ^^SPj 




...» 


; 8P„ 

7 


Si'^CP^ 

r’ 


— 


^ 1 ^CV 
1 


«;'-vp. 






7 


7*1 


— 




^ 2 ’^cVfj 








; 


•• 


: 


\ 


2 

“w' ^cv, 

r’ ' 


2 




Z 7 


7 ^*1 






V^CP^ 




2*' J 


®n’ ^8P|y 
7* ' J*) 




««',^BP„ 
J N 










; 




; 






.... 






— 


®n'^BPn., 










I 

1 ^ 


r*' 









where each block is of dimension R by R. There are N by W 
such blocks in (D-E)"^S. The following regularities in the 
block structure of (D-E)“^S are observed: 

1) Each block in the row of blocks has 8 / ^ as 



A. 



an operator where I is the index for . 



;th 



2) Each block in the j"*' ‘ column of blocks has either 



'BV 



. or S 



, as an operator where j is the index on wy • 



3) The upper-right and lower-left quadrants of 



(D-E ) are composed of blocks, each of which includes an 
operator of the form ^CPy- 

4) The upper left and lower right quadrants of 
(D-E)”^S are composed of blocks, each of which includes an 
operator of the form ^BV ;• 

^ _ j 

The matrix form of each block of (d-e) S is one of 



two kinds, B^Sj 0 p , or 

The form of is 



31 




B" ‘Sgp 



J 



( 31 ) 



where e^ and are values of equations (16c) and (16d) 
evaluated at | y .- | and co; is the weight in the quadra- 

ture set. 

The form of . is 



-32- 



U) o 







( 32 ) 



The fact that there are, at most, two terms in each 
element of (d-e) ^ S is attributed to the form of and 

^CVj each of which have, at most, two non- zero elements in 
each column. These properties give an ordered regularity 
to each block of (d-e)”^5. 

Equivalence of S/j Algorithm to Matrix Formulation 

As far as the author can determine, the S/j algorithm 
has not previously been described explicitly in the matrix 
form of the preceding section. The discovery that the Sj^ 
algorithm, described by equations (18) is equivalent to 
solving matrix equation (29) at each iterate step should be 
no surprise since equations (18) were derived directly 
from equations (17) . At each iterative step the Sj^ algo- 
rithm provides a recipe for inverting (D-e) and multiplying 
the result into 



and 2_. 
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This discovery, however, allows the Sfj algorithm to 
be classified as a method of successive approximation which 
provides a mathematical description of its convergence 
properties and error analysis. That is, equation (29) is 
precisely the same as (6). Consequently, properties 
previously discussed for the method of successive approxima- 
tions apply to the Sjy method if the conditions of Theorem 1 
are met. Before this is attempted, it is necessary to show 
that the norm of the iteration matrix can be calculated for 
the Sfj problem previously described. 

Calculation of | | (D-E)"^S| | 

The norm imposed upon the system is the infinity 
norm, defined for a vector by 

I lil L = . (33) 

k 

Obviously, this norm is conveniently chosen for pointwise 

5 3 

error analysis. Isaacson shows that for a matrix. A, 

I |A| loo = , (34) 

X. j 

that is, the maximum row sum of absolute values of the 
matrix elements. From this point, the above definitions 
will hold for the norms simply denoted by | | | | . 

The row sum of (D-Ej'^S, described by (30), (31), 

and (32) will be a row sum in a row of blocks which com- 
prise (D-E)"(s. Since each row of blocks has the same 
form, it is sufficient to sum a row of blocks and pick the 
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row from the block sum which gives the maximum row sum of 
absolute values for (D-E)'^S. 

Due to the ordering of the direction set, 



for j 
that 




( 35 ) 






,N/Z. 



A property of the direction set is 



W/2 

y W; • I . (36) 



Consequently a block row sum of (D-E)'^S can be collapsed 
by taking 



,-1o 



^CVfj '^CV^ 

-+ J il+ 2 

Z Z 



^ r - ^ ^ 



+ (8'.'S„^ +6'.'S^^ ) + . . . + (8,'So^ +B,'S(,p ) 






■ ^(W,*Wz*...*Wwl (Si'sBp* 82 'Scpl 



4 ' 1 



^IBl’Sgp.Bl’ScpI 



(37) 



where B.-^Sar) + B-^Si 



'X, ■’BP ^CV ~ matrix parts of equations 

(31) and (32). 



It is observed from (37) , (31) , and (32) that the sum 

of absolute values of the elements of each successive row 
include all elements of the preceding row plus a positive 
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quantity. This is true of all row sums except the last. 
The last row sum is equal to 



a 

r~ 



■i r 



i-i 



4-' 



+ 4 



l^-cl 



4 









( 38 ) 



The next to last row sum is equal to 



T 



The last row sum is the maximum if, 

R-2 





eR-2 


/ 










2 




R-3 


\] 




t- » 












el 








3 




X A 

'1 










4 


+ • • • + 







(39) 



R-2 



+ 3 



4-’ 



1 



or if 



+ 3 



R- / 



'/C 



(40) 



Even if A, a , and 1 | were adjusted so that < 0 , the 

quotient ^ I is less than unity because |e^| < d^ for all X. 
di 

Consequently, for systems with large R, the next to last 
row sum is the maximum. 

Now the block row sum must be picked from all possi- 
bilities. The index -c for which expressions (38) and (39) 
are maximum can be determined from equations (16c) and (16d) . 
The quantity 

m . Ji 






attains its smallest value for a given A and when |y| 
is the minimum in the quadrature set {yy}. Also 
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? 



2A 



is a maximum for the same | 



The 



- 2|P^I * Ao-^ ‘ 

latter factor dominates, however, in the norm determination. 
Therefore the norm is calculated by using the minimum 
absolute ordinate value in the set ^My)• Consequently 



(D-El’ 's I 




R- I 



3 


1 


S 7 I 


•J 7 


[i^i» 

|<3y 1 




» 






R-2 






2 


» 




z 




1 yf 




®^| 1® 








J 

37 




1 


d-r Id 




+ , . . + 


HZ 



R-2 






,R-2 






> ; 






otherwise 



(41) 



where ^ is the index for the minimum absolute ordinate 

value in the set {y>}. The latter value is appropriate 

1 

for large R. Observe that equality occurs when . 

For a given A,a'^,a'^, \ui\fnZn' ^ infinity- 

norm of the iteration matrix can be calculated by (41) . It 

is obvious that the norm is bounded since 1 < 1 and 4— 

is bounded for any finite dimensional system. 

Observe that, for a fixed A,a"^, and 

coefficients s, and are fixed. In this case the norm 
J i 

decreases linearly with a'^ . Consequently the norm of the 
iteration matrix is small for problems involving primarily 
absorbing media and large near unity for predominantly 
scattering media. This results in rapid convergence to an 
acceptable iterative error for absorbing media and slow 
convergence to the same level of error for predominantly 
scattering media. This behavior is predicted by equation 
(9) . It will be verified experimentally later. 
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Sufficient Conditions for Convergence of Sfj 

Equation (41) will be used to provide sufficient 
conditions to guarantee that ] ] (D-E)”^S| | < I in order that 
the algorithm can meet the conditions of Theorem 1. For 
large R, | | (D-E)~^S| | is evaluated by the second expression 
in (41), which can be rewritten 



(D-E)-^S| I < 






,R-2 



(42) 



where ^ 3 = 1 < 1 . 

For large R the terms ^ and 3^“^ 
could be neglected. But since 



(43) 

are insignificant and 




f 



imposing the condition 



a 



6 




< 1 



(44) 



will guarantee that | | (d-e)~^S| ] < 1 without seriously 
loosening the resulting bounds on A. 

From (16c), (16d) , and (43), expression (44) becomes 



or 



2\vi\ 



m^n 



Ao^ > 



Ao'* 



2 I y ; 



^ ‘ nu,n 



- Aa' 



(45) 
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This inequality provides the sufficiency conditions for 
convergence of the algorithm. 

Suppose the problem is defined so that ,o , and 
been selected and it is desired to impose a 
sufficient condition on the choice of A so that (45) is 
satisfied. Then the condition | | (d-e)'^S| | < I is satis- 
fied. 

Suppose further that Aa^ > so that < 0. 

Then (45) becomes 



or 



‘ — ■ 

o 

This imposes a limit on the negativity of the elements 
If had been chosen equal to Then 

CT ^ 

in (42) and -=— < / is a sufficient condition for 
I I (D-E)"^S| I < /. But this leads to the inequality 



(46) 



'-C' 



3=0 



^ ^ I I min 

2o^ - 



But ^ 




min 



so the inequality becomes — < — ; — or o'* < Conse- 

^ 2a^ - 

quently ||(D-E)"^S|| < 1 in this case if 



0 




< J . 



(47) 
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Finally, suppose Ac \ n\ln selected. This 

choice guarantees that ey and dy are positive for all j, 
hence (D-E)'^S and (D-E)"^ are positive operators. Expres- 
sion (45) is then 

^ 2Aa"^ + 2 1 Vi-4 1 iD-tn " 
or 

< ? . (48) 

Consequently, || (D-E)~^S|| < / for all physical problems 
with this choice of A. 

Properties of (D-E)"' 

The matrix (D"E) exists in the construction given 
by (21) and (22) . It is a linear operator due to its con- 
struction from the discretized linear equations (17) . It 
is a bounded operator since 

I I (d-e) II = dx + |e^| 

which is bounded for all finite values of cr^ , | | ^^^, 

and A . 

The matrix (D-E)”^ exists if cfe-t(D-E) is non-zero. 
From equations (21) and (22) 

de^(D-E) = (djd^- . 

and since dy > 0 for all j, dei(D-E) Jp 0. 
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All conditions of Theorem 1 are met by the 
algorithm if the sufficiency conditions listed in Table I 
are met by the problem parameters. If these conditions 
are met, | | (D-E)"^S| | < 1 and the Sj^ algorithm converges to 
, the exact solution of the discretized set of equations. 



TABLE I 

SUFFICIENCY CONDITIONS FOR CONVERGENCE OF Sj^ 



Possible Choice 
For A 


Value of 


Bound on A 
to Meet 
Sufficiency 
Condition 


Values of 

for Which 
Sfj Algorithm 
Converges 




Negative 

0 

Positive 


^ ^ 2 1 1 mZn 


0 < — < 1 


0-t 

A _ 




0 <^ < 1 
0 

0 < ^ < ] 


■ 

^ ^ ^ 1 1 nt'in 






at 



The above table can be summarized by stating that 

the sufficiency condition is met for any A provided that 

O'* < or iJjjLljZLiil, whichever is smaller. 

A 

The consequences of the algorithm meeting the 
sufficiency conditions for convergence are that expressions 
(1) through (9) can be used to study its convergence and 
iteration error properties. 
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Positivity 

Positivity of the solution of the Sfj equations for 
positive source and non-negative boundary conditions is 
desirable because it is a numerical approximation to a non- 
negative physical quantity. In practice the algorithm 
has produced negative solutions.® One-dimensional discrete 
ordinates codes such as ANISN^® and DTF-IV^^ have flux 
fix-up routines incorporated so that negative solution 
values are not allowed to occur. This strategy, although 
of practical merit, is not mathematically rigorous. The 
effect of these arbitrary fix-ups on the iteration error 
is not clear. Lathrop^® investigated various two-dimensional 
difference schemes and found that the quest for positivity 
is made at the expense of accuracy and efficiency of con- 

C C 

vergence. Bell points out that the Sfj difference equations 
may lead to negative solution values since they do not 
correspond to a positive operator. 

The preceding section revealed that sufficiency 

conditions for convergence of the S/^j algorithm can be met 

which allow the iteration matrix to possess negative matrix 

elements. But convergence to the exact solution cannot be 

guaranteed if A is chosen so that A > Conse- 

o 

quently, there is a limit on the negativity of some of the 
elements of (D-E)~^S if the theory is to be invoked to 
guarantee convergence. 

For computational economy reasons it is desirable 
to choose A as large as possible. For the problem 
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previously described, however, A < - I — | imposes an 
absolute upper bound on A , above which the sufficient condi- 
tion for convergence is not met. If this bound on A is 
observed, it is reasonable to assume that the approximate 
solution will be non-negative if -i is large enough 

even though initial iterates may have some negative com- 
ponents. This is because is expected to be a non- 
negative vector due to its approximation of ipj(x) which is 
non-negative by its physical interpretation and 
approaches as -i becomes large. This is true even though 

the theory of positive operators, which is used successfully 
in multi-group diffusion theory, is not applicable in 
this instance. Of course no guarantee exists that , the 
discrete approximation to ^j(k) is non-negative. Numerical 
experiments by the author confirm that it is for uniform 
and first collision sources when the sufficiency conditions 
are met. 

If A is chosen so that A < however, 

a 

(D-E)~^S and (d-e)"^ are non-negative matrix operators 
which guarantee a non-negative approximation for all Z. 



Necessary Conditions 

The previously discussed conditions are not necessary 
for the Sjy iterative process to terminate based upon 
fractional differences reaching an arbitrarily small value. 
Indeed termination of the Sj^ algorithm iterative process 
has been observed by this investigator for choices of 
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^ 1. ^^ I 'Pj' ! } . for which ||(d-E)"^S|| >> 1. For some sources, 

e.g. uniform source, convergence to a positive solution 
vector occurred. For other sources, e.g. first collision 
source, the solution vector contained negative components. 

It is doubtful whether this process can accurately be called 
convergence since the nature of the resulting solution is 
not clearly defined. 

Even if A is chosen so that A > ^ I . ^A , l , in which 

at 

case all ey are negative, (d~E)'^ 5 contains even and odd 
powers of Consequently, the iteration matrix is not 

wholly negative. In this instance the S/j algorithm may or 
may not converge to a non-positive solution depending upon 
the source distribution. For these cases, where some 
elements of the solution are negative, doubt exists as to 
the nature of the solution. 

It is not clear that necessary conditions on A exist 
for convergence of the algorithm. It is clear, however, 
that a solution exists if the sufficiency conditions are 
met and that these conditions allow a limited amount of 
negativity for the iteration matrix. For the numerical 
problems previously cited, these solutions are positive. 




m 
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Effect of A on I I ( D- E) ' I I 

Figure 2.2 reveals the 
effect of A on | | (D-E)"^S| | 
for the following problem 
parameters : 

slab width = 10.0 cm 

^ « -1 
0 = 1 . 0 cm 

= 0.9 cm~^ 

= 0.23862 

As R (the number of spatial 
increments) decreases, A 
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A, c»-v. 



becomes larger and | | (D-E)~^S| | increases. In the region 
where A meets the sufficiency conditions, | | (D-E)”^S| | < 1 



and convergence of S/^ to the exact matrix solution is assured. 



A small part of this region allows negative elements in 



(D-E) ■ h. 

In the region where | | (D-E) | > 1, termination of 

the iterative process may occur. Solution vectors may be 
positive or partially negative depending upon the source 
distribution. This region should be avoided for practical 
problems unless other evidence is available that the solu- 
tion has validity. 



Effect of Quadrature Set on | | (D-E)"^S| | 

Expression (9) reveals that the iterative error is 
less tightly bounded as | | (D-E)”^S| | increases towards 
unity. This means that more iterations are required to 
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reach a specified accuracy. Equation (41) shows that 
||(D-K)"^S|| increases for fixed , as the ordinate 



|M(|m(>( decreases. This is because ^ grows larger and 
|e^-|/d^- diminishes as A increases but predominates in 

(41) so I I (D-E)“^S| I increases. All quadrature sets have 
the property that approaches zero as the order W of 

is increased. Consequently, the effect of choosing a 
higher order is to increase | | ( D-E) “ | | , hence requiring 

a larger number of iterations to convergence. This loss of 
efficiency might be observed in going from to where 
the spread in \)^i\min values is significant but is of no 
practical importance for higher order Sf^ where \ Vi] min 
values are quite close. 



Convergence Properties of Sf^ Algorithm 

The following discussion applies to the Sf^ algorithm 
when A < A I . Then e . and d - are positive for all j 

0 -t j j 

and (D-E)”^S and (D-E)~^ are non-negative. The sufficient 
condition for convergence is met for all physically 
realizable problems. The initial guess, - 0, is a con- 

venient reference so that comparison of convergence 
efficiency can be made with another method later. 

The algorithm, under these conditions, converges 
in such a way that each vector component approaches its 
exact value in a monotone non-negative manner. This can be 
ascertained from applying equation (7) to the difference 
defined by 
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( 49 ) 



Then 



■c 

1 

m-o 



■i- 1 



(d-e)"^S]'”(D-E)‘^£ - ^[(D-E)"^S]"'(D-E)"^£ 



m~o 



= [ (D-E) - ^S]^(D-E) - 



(50) 



which is positive for ^ > 0. Consequently, some positive 
quantity is added to each component of to produce the 

respective component of 

From the definition of given by (8) and given 

by (4), 



■i- 1 



= 5! t (D-E) ' ^S] (d-E) “ ^2.- [ (D-E) ■ ^5] *" (d-E) ” 



m= 0 



m= 0 



= [(d-e)"^S]^{I+(D-E)'^5+[(D-E)-^S]^+---}(D-E)‘^£ (51) 



which is positive for ^ > 0. But 









(52) 



Hence 



jl-cMI , ^Ul 



(53) 



due to positivity of and See Figure 2.3 for a 

convergence diagram. 
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The largest absolute error in each successive iterate 
is monotone decreasing. From equations (51) and (4) , 

= [ (D-E) ■ . (54) 

Consequently , 

= [ (D-E) - = [ (D-E)- ^S] . (55) 

This means that 

_ {(q_£) for all k which 

k - ri 

n (-C ) 

does not imply that is some fixed fraction of 

for all k even though (52) ensures that ^ ^ for 

all k. That is, it cannot be shown that ^ ^ < 3 

for all ■i,k where 3 is a positive constant less than unity. 

The most that can be inferred from (55) is 

< IKD-EI-'sll.lle'-^'ll • (56) 

This means that the maximum error at each iterative step 
diminishes if | | (D-E)”^S| | < 1 . This condition is met for 
convergence . 

Similarly it can be shown from equation (50) that 

I I I < II (D-E)'^5| 1 -I 1^*^’ I I (57) 

using the same arguments and with the same consequences 
discussed above for the error. 




p' 
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It is not of much practical value to know that the 
error at each iterative step is less than the maximum error 
at the previous step, it would be of more practical value 
to know that the error is bounded by some fraction of the 
difference at each step because the differences can be 
calculated. Equation (51) is 

= [ ( D-E) ■ (D-E) ■ ^S+ [ (D-E) ■ ^S] ^ + * • • } (D-E) ■ 



= {I+(D-E)"^S+[ (D-E) " ^S] ^ + * ♦ • } [ (D-E) ■ ^S]'^(D-E) ■ 

because [|D-E)~^S]'^ commutes with each term in the braces. 
So 

= {I+(D-E)"^S+[(D-E)'^S]^+*-*} 



due to (50) . Consequently, 



e 



U) 

k 



< 




II ^ 



iii''^*’'ii 

I- 1 I (D-E)'^S| I 



(58) 



for all k. This does not provide enough information to 
assure that is some fixed fraction of for all fe 

but it is enough to provide a practical link between the 
error and difference which can be exploited to derive a 
useful convergence criterion. 
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Convergence of the Sf^j 
algorithm is pictorially 
described in Figure 2.3 for 
each vector component. The 
iterative solution proceeds 
progressively towards the 
exact solution in a monotone 
positive way governed by 
properties (53) , (56) , (57) , 

and (58) . 
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Conventional Pointwise Convergence Criterion 

The pointwise convergence criterion based on 
fractional differences. 



,(^+n ,u) 

h -♦fe 



W) I 'C ) 
k 



< 



e 



(59) 



where e is chosen arbitrarily small, is used universally in 
practice for shielding problems. It is especially useful in 
deep penetration problems where has elements of very 
small magnitude since it weights these regions more heavily. 
That is, iteration proceeds until the deep fluxes, which 
are the last to meet the convergence criterion, are con- 
verged. Its utility is based on the practical idea that, 
as the differences become smaller with successive iterations, 
the fractional amount added to each solution component is 
smaller than the desired accuracy after some Consequently, 
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further iteration will not add significantly to the solu- 
tion. 

It must be ascertained, however, whether this con- 
vergence criterion, if met, implies that the pointwise 
fractional error is sufficiently bounded. After all, it is 
desired to iterate until assured that the pointwise iteration 
fractional error is arbitrarily small, not just until the 
pointwise fractional differences are bounded. 

Suppose iteration proceeds until condition (59) is 
met. Since 

this implies that 



< e ^U) 



but = (D-E) so the above implies 



(D-E) ■ < e 



or 

< e + (D-E)'^S . 



But by definition 



[ (D-El'^Se^-^hfe ^ I I (D-E)- I I 



( 60 ) 



< I I (D-E)'^S| I • I I I . 



-51- 



Consequently (60) implies that 

e^l < e , I I (d-e|-’s| l-l | | 

for all k. But from the definition 

^(•^1 = the above inequality becomes 

4 ^' < e ♦ II |D-E|-'s||.||e'-"'|| 

or 

< eifj ♦ I I (D-E|-'s| | -| ll''"’l I 

for all fe. But if this inequality holds for all k, it 
must hold for j, the value of k which corresponds to the 
maximum absolute value of - Hence 

I I I I * [r+e-l I (D-E)-^sl l] < e ^ j* 

or 



J+e-| I (D-E)'^S| I 

where j is the index for which attains its maximum. 

Observe that (61) does not imply that the fractional error 
is bounded for all components of \p * - However, if it is 
arbitrarily assumed that 



ip[^* > ipj* 



( 62 ) 



for all k, expression (61) implies that 
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■7*“ 



J + e- I I (D-E)-?S| I 



or that 



l-i) 

< f (63) 

W )+e- I I (D-E)-?S| I 



for all fe. But even under the dubious assumption (62) , this 
does not imply a satisfactory bound on the fractional error. 
Figure 2.4 displays a plot of the bound versus | | (D-E)”^S| |. 
Observe that, as ||(D-E)"^S|| approaches unity, the bound 
grows to 1. That is, up to 
100 percent fractional error 
is implied by invoking (59) as 
a convergence criterion. This 
criterion is useful and 
practical for | | (D-E)”^S| | < 

0.9 but for larger values it is 
a poor choice because the 
resulting bound on the 
fractional error is too large. 

An Improved Pointwise Convergence Criterion 

Implicit in the concept of a practical convergence 
criterion is that iteration is terminated only after the 
fractional error is arbitrarily small. Since the exact 
solution is not known, it is necessary to iterate on 
fractional differences, which are known. | | (D-E)“^S| | can 
be calculated for the stated Sj^ problem so it is possible 



FIGURE 2.4 

BOUND ON FRACTIONAL ERROR 
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to establish a convergence criterion based on fractional 
differences which implies that the fractional errors are 
arbitrarily small. 

Suppose iteration proceeds until 



6 n I 



ii^U. 

k 



yEc I (D-E)'^S| f 



(64) 



for all k where e is some arbitrarily small positive number. 
This implies that 



(^+/) 



e ijj 



k 



J- I I (D-E ) - 's 



j-e for all k. 



But repeating (58) 



U) 



U^l) 



T-l I (D-E)-^Sl 



so the above 



inequality implies that 



1] 



£ \U 



for all k. 



But by definition , so the above inequality 

implies that 

^(x.) < ^ 
or 



e 



(-t) 



< e ip* 



or finally 



IT 



U) 






< G 



(65) 
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for all fe. Consequently, upon invoking convergence 
criterion (64) , the convergence properties of the Sfj algo- 
rithm assure that the fractional iteration error is 
arbitrarily small. This is precisely the goal for a machine 
convergence criterion. 



as the convergence criterion but recognition of its inability 
to guarantee a satisfactory bound on the fractional error 
for problems involving strongly scattering media is impera- 



0.9, the conventional convergence criterion is acceptable. 
Table II shows that the effect of using the new convergence 



For the Sfj problem discussed here, equation (41) 
allows the exact computation of | | |D-E)‘^S| | and e is 
arbitrarily chosen. Of course for problems in which 
||(D-E)"^S|| is not obtainable, it is necessary to use (59) 



tive. 



FIGURE 2.5 



Figure 2.5 reveals that 



BEHAVIOR OP IMPROVED CONVERGENCE 
CRITERION 



the convergence constant. 



j^/ - I I { D- E) ” I ij becomes 




c 

i-e 



vanishingly small as 
||(d-E)"^S|| approaches unity. 
For predominantly scattering 
media ||(D-E)~^S|| is quite 



U/ I 

close to unity. It is neces- 



sary to use the improved con- 



||(D-E)-‘S|1 



1 



vergence criterion (64) for 



these problems. For weakly scattering media, i.e. a'^/a^ 



< 



criterion is small for 




< 0.9 but overpowering 



for I I (D-E)"^S| I > 0.9. 




% 

f ■ 



5 '. ■ 

■f 
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TABLE II 

IMPROVED CONVERGENCE CRITERION 
FOR e = 10"4 



1 1 (d-E) ■ I I 


[l-||(D-E|-IS||] 


0.5 


0.5 X 10“"^ 


0.9 


10"5 


0.99 


10"^ 


0.999 


10"'^ 



The improved convergence criterion imposes a much 
more severe condition for convergence than does the con- 
ventional convergence criterion, especially for predominantly 
scattering media. For these problems, convergence is slow 
even using the conventional criterion. Imposing the 
improved convergence criterion makes a bad situation worse. 

It is necessary, however, if assurance is required that the 
fractional error is arbitrarily small. 

Computer Experiment 

A sample Sjy problem, described by equations (18) , 
was run on a 360/67 IBM computer with the following param- 

Slab width = 10.0 cm 
R = 30 spatial intervals 
= 1.0 cm“^ 

= varied 
e = 10“^ 



eters: 
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N = 6 ordinates using Gauss quadrature set^^ 

Uniform source, Qylx) = 1.0 for all /. 

It is worthy of note that any of the more recently derived 

2 3 62 

mechanical quadrature sets, ' could have been used but 
the Gauss set is satisfactory to demonstrate convergence 
properties of the algorithm. This set will be used 
throughout. 

A copy of the computer program and samples of the 
computer output are listed in Appendix A. Table III lists 
the results. The number of iterations to convergence are 
tabulated for various a'^/a^ ratios using each of three 
different convergence criterion. 

Observe that, regardless of the convergence criterion 
used, problems involving predominantly absorbing media 
required fewer iterations to converge than did those for 
scattering media. The number of iterations to convergence 
is a monotone increasing function of 

Comparing columns two and three shows the effect of 
converging on the maximum difference instead of the point- 
wise difference. As expected, the maximum difference 
criterion is more severe since it takes more iterations to 
meet it than the pointwise difference criterion for each 
value of This is expected because 









UT 



U+ 1 ) 
fe 



for all fe by definition. 



NUMBER OF ITERATIONS TO CONVERGENCE 
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Comparison of columns two and four shows the effect 
of the improved convergence criterion which guarantees 
that the fractional iteration error is less than e. Its 
effect includes the effect discussed in the preceding para- 
graph and has the additional effect of | | { D-E) ” | | . Both 
effects increase the number of iterations to convergence. 
The effect of |1(D-E)“^S|| becomes more dominant as o'* / 
increases above 0.9. This is due to the direct linear 
relation has on | | (d-e) | . 

A comparison of solutions, one of which met con- 
vergence criterion (59) , the other meeting the improved 
convergence criterion (64) , was made for the values of 
0 - 6 / 0 -^ listed in Table IV. 



TABLE IV 

COMPARISON OF SOLUTION VALUES 

e = 10"4 



/a^ 


^ ; 1 

Agreement of Conventional Convergence Criterion So- 
lution with Improved Convergence Criterion Solution 


0.5 


Agreement in 4th significant figure 


0.9 


All agree to 2nd significant figure, many to 3rd 


0.99 


All agree to 2nd significant figure, some to 3rd 


0.999 


All agree to 2nd significant figure, some to 3rd 


0.9999 


All agree to 2nd significant figure, few to 3rd 



For this problem ||(D-E)"^S|| = a^/a^. 
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The above table reveals that agreement exists only 
to the second significant figure for o^/o^ > 0.9. From 
previous discussion it was established that a more accurate 
solution results when a larger number of iterations are 
made. Consequently the improved convergence criterion solu- 
tion must be closer to the exact solution than the con- 
ventional convergence criterion solution because more 
iterations were performed for the former. The above table 
shows that the conventional convergence criterion solution 
agreed to only the second significant figure with the more 
accurate improved convergence criterion solution. This 
illustrates that the expected accuracy of e = 10“'* was not 
achieved for large values of o'^lo^ using the conventional 
convergence criterion. 

Tables III and IV provide vivid illustrations of the 
necessity for using the improved convergence criterion for 
problems with I1(d-e)”^S|| > 0.9 if an arbitrarily small 
fractional iterative error is required. 

The overall effect of using the improved convergence 
criterion is insignificant for problems involving strongly 
absorbing media for which | | (d-e)~^S| | is small. The effect 
is marked, however, for predominantly scattering systems 
and results in additional computational work required to 
assure arbitrarily small fractional iterative errors. Con- 
sequently any scheme which accelerates convergence of the 
algorithm by decreasing | | (D-E)”*S| | would be eminently 
useful when invoking the new convergence criterion, especial- 
ly for strong scattering media problems. 
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CHAPTER III 



ACCELERATION OF Sfj CONVERGENCE BY SPATIAL TRANSFORMS 
General Discussion 

The previous chapter discussed how the Sfj algorithm 
converges more rapidly for absorbing than for scattering 
media. As decreases, holding constant, the number of 
iterations to convergence decreases regardless of the choice 
of convergence criterion. Based on this physical argument, 
it is expected that a mathematical device which adds 
absorption to the problem will accelerate convergence. From 
a mathematical point of view, an operation on the 
equations which decreases the norm of the iteration matrix 
will improve the efficiency of the algorithm. A new method 
is presented here which incorporates both ideas. 

General Spatial Transform 

A transform of the type 

'|»U) = <}>{x) (1) 

has been used to provide an effective importance sampling 
device in the solution of the transport equation by Monte 
Carlo methods. It biases the sampling distribution 
towards more important directions, thereby reducing the 
variance in Monte Carlo methods. This leads to improved 
statistical accuracy for a given calculational effort. 
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Transforms of this type have been used by the author 
to accelerate convergence of the algorithm for certain 
problems. Transforms which have proven successful are: 



'Pj 


(x) 




(x) 


ax 


(2) 




(x) 


- 


(x) 




(3) 














*J 


(x) 


- 


(x) 


r x+b 1 

( x+b ] J 


(4) 



where a is an arbitrarily chosen acceleration parameter. W 
and b are parameters chosen to restrict the range of the 
arguments. These transforms are applied to the angular 
discretized transport equations. The most effective trans- 
form of those listed is (2) . 

The Half-Range Problem 

The spatial transform accelerates convergence of the 
algorithm for the deep penetration infinite slab problem 
under the following assxamptions : 

i) steady state 

ii) homogeneous medium 

iii) angular flux has no azimuthal dependence 

iv) high energy particles are scattered only in forward 
directions 

v) forward scattering is equally distributed in 
direction. 

A consequence of these assumptions is that the trans- 
port equation has the form 
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.Slip. 

^ 9x 



+ a 



^ 1\1 I } ctu ' + <?(x,y) (5) 



valid for 



0 < y i J . 



The anisotropic forward scattering assumption is valid for 
high energy neutrons and gammas. The assumption that the 
forward scattered particles are equally distributed in 
angle is a simplification of the usual Legendre Polynomial 
approximation and is made strictly for convenience. The 
more accurate description of the forward scattering aniso- 
tropy could be treated by the method. 

The scattering distribution function, under assumption 
v) is normalized by recognizing that ^(y'->-y) = aon6i and 

1 

J ( coyut) dy = / . (6) 

0 

This simply states that the probability that a scattered 
neutron emerge into some allowed direction is unity. A con- 
sequence of this normalization is that 

= I . (7) 

Angular approximation is accomplished by discretizing 
equation (5) in the W/2 forward directions of the previously 
described quadrature set. This gives the approximation 

d^\>.[x) . W/2 

— dx — * 



(8) 






I 
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with j = Quadrature weights obey 



W/2 

I “n = ^ 

n= 1 



( 9 ) 



and all y ; > 0. The angular quadrature set is the positive 

j 

half of the one used for the full-range problem. 

Spatial discretization is accomplished in the same 
manner described in the previous chapter to give 



J^k+1 ,j 



e .ih , 
j 



4 W/2 
n= I 






(' 1^1 . ■'’'P , ) + Q 

H ^ n fe, K ^ 






(10a) 



for y = and k » 7, 2, •••,!? 

where 

2y,- + Aa^ 
d/ = _i 



and 



2y , - Aa^ 




(lOb) 



(10c) 



and R is the number of equi-spaced spatial increments, A. 

The algorithm, as before, solves this system of 
equations iteratively as described by 



i 



(-t+n 






oj 

2 



J n \ fe+ 7 , n fe, d/ 



fe+7/2,y 



( 11 ) 



where ^ is an iteration index. This set of equations is 
divided by dy to give the computer-solved algorithm 
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iUU 



^kj *dy 1 ^„f, ^k.n} 



^ "^fe+?/2,y 



( 12 ) 



Numerical results for this algorithm with the first 
collision source 



. = e 






and vacuum boundary conditions 



j/y . 



(13) 



^ 



(14) 



for all y, will be displayed and discussed later. 

Iterative Matrix Formulation 

The Sfj algorithm (12) , (14) with source (13) has the 

matrix description 



(D-E )i^ ^ ^ ^ 



(15) 



derived from equations (11) . 

If the equations are arranged so that the solution 
vector has the order 






4' 






3,1 






R+1 ,1 

2,2 



3,2 
R+ / , 2 



'P 



R+?,W/2 . 



( 16 ) 



- 65 - 



then 



"^ 3 / 2 ,; 



1 = 




^ 512,2 
"^ 5 / 2, 2 




( 17 ) 



The matrix (D-E) has the block form 




( 18 ) 



where each block has the lower bi-diagonal form 
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- e - 



d. 




- e . 



d. 









(19) 



The index Z is that for the ordinate Uy. Each block is an 
R by R submatrix. 

The scattering matrix has the form 



^BVj 




• • • • 

• • 4 • 




^BVI 


^BV^ 




• 




• 

• 

• 


^BVj 


^BV^ 


• • • « 





( 20 ) 



where each block has the lower bi-diagonal form. 







4 



i 



(21) 
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and 



6 . 
i 



6 




( 22 ) 



Observe that this matrix formulation is the upper left 
quadrant of blocks of the full-range matrix formulation 
and only the elements Ay are different. 

Due to the block diagonal structure of (D-E), and 
the bi-diagonal form of each block, (D-E)“^ is simply 





= (D-E)'^ (23) 




where 



68 



l^th 






ex 

2 


1 

dx 

eX 


dX 

• 


2 


dx 


J-1 

4 


CD 

1 


• 

• 


4-2 


# 

A-3 


dH-^ 




yL 


yC 




4'^ 




yC. 




4'^ 


yC 


yL 








(24) 




and /t is the index for y-. Each block is of dimension R by 
R. (D-E)”^ can be applied in its general form to equation 
(15) to produce 



= (D-e) ■ + (D-e)‘^£ • 



(25) 



In this equation 



(D-E) 











r“N 

“2 ^BPj 


^2 




• # 

• • 

• • 






^N/2^BVj 


R-^ S 

“W/2^8P2 







( 26 ) 
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Observe that this matrix is the upper-left quadrant of 
blocks of the full-range iteration matrix of the previous 
chapter. 

Each block has the form 




Calculation of the Norm 

The infinity-norm of (d-e)”^S, as previously defined, 
is easily calculated. The block row sum is 



^Bt) ^BD ~ +•••+(») )(B.^Sp„) 

t aVj ^ bv^ A, 2 1 2 NI2 -t- SP' 



rf^ - / 

= ?- B. S 



BP 



(28) 



due to (9) . Here Sg^ has the form of (21) but with unity 

a-* 

in each non- zero matrix element. Actually (28) is ^ 



times 
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thc matrix part of equation (27) . Consequently 



(D-E)"^S| 



R- ? 



R-2 



'■L 



+ 21 



iiL+ifit +■ 



e; 






L d 


















)J 



(29) 



where as before, -c is an index for the minimum absolute 
ordinate value in the set {yy}. If one defines a 6 



0 




3 < I 



(30) 



equation (29) becomes 



I I (D-E)-^S| I 




(31) 



Sufficiency conditions for | | (D-E)"^S| | < J are derived by 
requiring that 



^ 1-6 



? 



or that 

2|y .J + Aa^ > 2Aa^ 



(32) 



This is the same inequality which gave the sufficiency con- 
ditions in Chapter II for convergence of the full-range 
algorithm to the exact solution These conditions are 

the same for both problems and will not be restated here. 
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Exponential Spatial Transform 
The transform 



=4>y(x)e 



ax 



(33) 



where a is a constant greater than zero, applied to the 
angularly discretized equations ( 8) gives 



d<p-(x) W/2 

ly — ^ + apy)(l)y{x) = 



e (34) 



j = Observe that the transform has added 

absorption to the problem in amount ayy to produce a new 
"effective" total cross section, {o^ + ajjy) as well as 
changed the source . 

Spatially discretizing these equations in the same 
way as previously described for the unmodified equations 
gives 



J 2 k+l,n ^k,n^ 



+ 1 

r (35) 



k = and j = ] ,1 , ' ' ' ,W j 1 , The algorithm in machine 

form is 









(36) 



where 



ay. 

- a , ^ 



2y; + Ao + Aay - 

J J 





(37a) 



and 
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e. 

J 





- AaM . 

i 

2A 



(37b) 



Matrix Formulation 

The (|)-domain algorithm has precisely the same matrix 
form as the unmodified equations. That is, equations 
(11) through (27) describe the (f>-domain algorithm if ey and 
dy are replaced by g.y and dj and if source and solution 
vectors are changed accordingly. 

The infinity norm of the transform domain iteration 
matrix is 



I I (P-E)‘^S| I 




2 









(38) 



where i. is the same index as before. Defining 




The above equation becomes 



I I (P-E)" 'S| 





(39) 



(40) 



Properties of ||(P-E)"^S|| 

The infinity norm of the transform domain iteration 
matrix is a monotone decreasing function of the acceleration 
parameter, a. For small a, the slope of | | (P-E)”^S| | with 
respect to a is steep but 'decreases to near zero at a = a* 
where 
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a 



* = 






- Ao' 



A ly 



-c ' m-cn 



(41) 



which is the value of a when 2.^ is zero. The norm value at 
a* is 



I I IP-E 




a=a* 



Aa^ 

TWJl 



m^n 



(42) 



From this point the norm slowly approaches zero as a 
approaches «>. This behavior of | | (P-E) ^S| | with a is 

verified by analysis of equation (40) in Appendix B. 



Optimum Value of a 

Appendix B reveals the behavior of | | (P-E) | 

for a > a* to be 



I I (P-E)"'S| I = 



Aa^ 

FTiT" 



) - 



2 (Ae) 



R-r 



^ m^n 



(4 + Ae 1 



R 



( 43 ) 



where 



and 



a = a* + e 



(44) 



0 < e < °° . (45) 

In the region where Ae < 1, | | (P-E)'^S| | essentially takes 

its value at a*. Consequently very little norm reduction 
occurs in this region. This defines a practical "optimum" 
value for a, that is, a = a*. Figure 3.1 illustrates this 
"optimum" value for a particular problem. Observe that the 
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norm is not further reduced for values of a > 2.00, at 
least to five significant figures. For this problem, 

(X* = 1.80924. Hence an Sf^ problem of the type discussed 
here has a practical "optimum" acceleration parameter which 
can quickly be calculated from equation (41). Use of this 
value will result in the largest practical norm reduction 
with accompanying improved iteration efficiency. 

Range of Utility of Transform Method 

From the previous discussion, it is obvious that 
I I (t>-E)"^S| I has reached a practical minimum value when 
is zero. Consequently the transform method can signifi- 
cantly decrease the norm of the unmodified Sj^ iteration 
matrix only when (D-E)"^S is non-negative, that is, when 

^ min . (46) 

Otherwise is negative, which guarantees that is nega- 
tive for all a >. 0. In this case, the norm has already 
reached a value which decreases very slowly with increased 

a . 

The transform method has the greatest latitude for 
norm reduction, if (46) is met, when A is chosen very small. 
Equation (41) reveals that a* is large when A is chosen 
small. A large a* decreases the norm of | | (P-E)”^S| | more 
than a small a*. Table V shows this property for a specific 
problem. 



I 

■#g 
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TABLE V 

a* FOR VARIOUS A 

Problem parameters: 
slab width = 10.0 cm. 
= 1.0 cm“l. 

= 0.9999 cm"^. 

(y.) . = 0.23862 

I I (D-E) ■ I I = 0.9999. 

Table V shows that small A allows large values of a 
to be used resulting in vast reductions in | | |d-E)”^5| | 
and improved convergence efficiency. The choice of a, how- 
ever is based primarily on the estimation of tolerable 
discretization error for a given problem. Consequently 
the transform method's utility depends upon factors other 
than those which most improve iterative efficiency. Once a 
A is chosen, however,' a* is easily determined. Its use 
provides practical maximum acceleration of convergence for 
that problem. 

Computer Experiment 

The half-range problem algorithm previously described 
in this chapter was run on an IBM 360/67 computer with the 
following problem paraimeters: 

Slab width = 10.0 cm 
R = 30 spatial intervals 
A = 0.3333 cm 
= 1.0 cm~^ 



R 


A 


a* 


1 1 (t>-E)-^S| 1 at a* 


30 


0.333 


1. 81 


0.69839 


50 


0.200 


5.81 


0.41903 


100 


0.100 


15.81 


0.20952 
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= 0.9999 cm“l 

Six ordinate Gauss quadrature set 
First collision source of unit strength 
Fractional iterative error less than e = 10“^. 

The improved convergence criterion, described in 
Chapter II, was used. A program listing and sample output 
for unmodified are included in Appendix C. The transform 
method program and sample output are included in Appendix D. 

Table VI compares the number of iterations to con- 
vergence for various values of the acceleration parameter, 
a. 



TABLE VI 

TRANSFORM METHOD RESULTS 



Alpha 


Infinity 

Norm 


Number of 
Iterations 
to Converge 


Percent 
Decrease 
in Number of 
Iterations 


Maximum Percent- 
age Difference 
in Unmodified to 
Transform 
Solutions 


0 


0.99990 


59 


— 


0 


(N 

• 

O 


0.95435 


48 


18.7 


0.36 


0.5 


0.89332 


46 


22.0 


1.1 


1.0 


0.80727 


45 


23.8 


7.1 


1.809* 


0.69839 


45 


23.8 


38.9 


2.0 


0.69839 


45 


23.8 


49.4 



Observe that a = 0 is the unmodified Sjy/ result. The 
symbol * signifies the optimum a. For R = 30, the condition 
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‘ ' ^‘■a 'mJ n 

A < ^±1. is met and the unmodified Sj^ iteration matrix is 

non-negative. 

Table VI reveals that the transform method accelerates 
convergence of the Sfj algorithm for the half-range problem. 
The acceleration parameter, a, decreases ey and increases 
dy so that each non- zero element of the iteration matrix 
is decreased in value (see equations (26) and (27)). This 
produces a reduced norm and concomitant acceleration of 
convergence. 

Values of a greater than a* used here do not further 
reduce the norm or number of iterations to convergence. 

The dependence of convergence efficiency on the norm is not 
linear. Acceptable improvement of convergence efficiency 
is attained for values of a well below a*. 

Calculational Effort 

The calculational effort of a machine algorithm is 
based upon the total ntamber of machine operations performed. 
Since additions and subtractions require only an infinitesi- 
mal amount of computer time compared to multiplications 
and divisions, only the latter are tabulated. Table VII 
lists the number of operations in each component part of the 
unmodified Sjy and transform method algorithms . Here n is 
equal to (R*W/2), -i is the number of iterations to convergence 
of the unmodified algorithm and J is the number of itera- 
tions required for transform algorithm convergence. 



-79- 



TABLE VII 

COMPUTATIONAL EFFORT COMPARISON 



Quantity 


Unmodified Sjy 
Algorithm 


Transform Method 
Algorithm 


Mesh sweep 


-t* 1 2n) 


y- [2n) 


Scattering source 


-c[W/2 + 2] 


y [W/2+2] 


Coefficients e, d 


9IW/2) 


131N/2) 


Transform source 


0 


n 


Transform solution 


0 


n 


Total 


l[2n+^+2]+9 ij) 


y[2n+^+2]+2n+/3(|) 



The transform method reduces computational work if. 



j < I 



4{3N-1] . 

wnwTj 



(47) 



For most problems R >> N so the transform method has utility 

even if only a few iterations are saved in its employment. 

53 

Isaacson shows that Gaussian elimination requires 



n[n^~ 1 ) 

3 



+ n 



(48) 



operations to solve a matrix system corresponding to the 
exact solution, Comparison of the Sf^ iterative technique 

shows a significant economy for large systems, n large, if ^ 
is not too large. Even so, I would have to be of the order 
before the computational effort of the unmodified Sf^ 
algorithm reached that of Gauss elimination. 
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Error in Transform Method Solution 



The new convergence criterion 



max 

y 




[l-l I 1P-E|-'S| l] 



(49) 



which guarantees that 




(50) 



for all k, was used in the computer algorithm. Using the 
transform (33) / inequality (50) implies that 



solution. Consequently the convergence criterion (49) and 
convergence properties of the method guarantee that the 
transform method i|;-domain solution has fractional iteration 
errors less than e. This same criterion was met for the 
unmodified Sjy algorithm. The origin of the discrepancy in 
the (|;-domain solutions for transform method and unmodified 







or that 




< e 



(51) 






for all k, where is the (/j-domain transform method final 
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Sf^ solutions, exhibited in column 5 of Table VI, is there- 
fore not attributable to iteration error. 

The origin of the discrepancy mentioned above is 
discretization error. The (j>-domain exact solution, ^*, has 
much greater variation over the slab width than does the 
unmodified Sjy exact solution, ij^*. Figures 3.2 and 3.3 
demonstrate this for a particular direction. Since the 
spatial discretization of the derivative is accurate to 
order in both domains, it is a poorer approximation in 
the <j)-domain, where the solution has large curvature magni- 
tude, than it is in the unmodified Sfj i|^-domain where the 
solution has a relatively small variation over the slab 
width. Consequently larger discretization errors occur in 
the transform method (J)-domain solution when large values of 
a are used. At ot*, varied seven orders of magnitude 

over the slab width whereas the unmodified Sj^ i|)-domain 
solution varied by less than a factor of two. This accounts 
for large discretization errors in the i|>-domain transform 
method solution when is near a*. Table VI reveals, how- 
ever, that as a is decreased from a*, the discretization 
error diminishes so that the i|)-domain solutions for both 
algorithms agree to within less than 0.4% at a = 0.2. For 
this a, the c(>-domain solution varied less than a factor of 
five over the slab width, producing smaller discretization 
error than did larger oi. Here convergence was accelerated 
to 78.5% of the maximiim acceleration obtained at a*. The 
value of a = 0.2 was a good practical choice for this 
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FIGURE 3.2 

SAMPLE UNMODIFIED Sfj >1' -DOMAIN SOLUTION 




X.CM. 




> 




NORMAL I ZED ANGULAR FLUX 
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FIGURE 3.3 

S7VMPLE TRANSFORM METHOD 4»-DOMAIN SOLUTION 




8 .. 



8 .. 

8 .. 




Parameters : 
a = a* = 1.80924 
R = 30 
A = 0.333 
= 0.9999 
= 1.0000 



U.OO 6.00 8.00 10.00 

X,CM. 



2.00 
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problem. For other problems, the user must make a choice 
of a based on practical experience and a guess at the 
expected solution. 

Influence of A 

Previous discussion connected with Table V revealed 
that, in theory, a smaller A allows the choice of a larger 
a* with accompanying larger decrease in | | (P- E ) ” | | and 
number of iterations to convergence. In practice, however, 
such large values of ^ are not useful because of the 
extremely large resulting discretization errors in the 
(j)-domain solution. Another consideration also limits the 

practical usefulness of large cx. For thick slabs, the trans- 

- QtX 

formed source, which has a factor e , can become extremely 
small. Since all computers have a limit on the smallest 
number it can handle, 16~^^ for the IBM 360/67, it is easy 
to exceed this limit for large a. 

As a consequence of the above factors, it is best to 
use values of a much smaller than a*, for a problem with any 
A for which the transform method is applicable. Computer 
experiments conducted by this investigator have verified 
that the choice of a= 0.2 produces approximately the same 
acceleration of convergence and accuracy listed in Table VI 
for the problem described previously but with the much 
smaller values of A resulting from R =* 50 and R = 100. 




I 
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Generalization of Transform Method 

The transform method immediately generalizes for the 
half-range problem treated here to include heterogeneous 
slabs. In this case the cross sections are functions of 
position. The relatively simple expressions for ||(D-E)"^S|| 
and I I | shown in this chapter do not occur but the 

matrix block form prevails and these norms can be calculated. 

Transforms such as equations (3) and (4) place space 
dependent absorption into the problem. When using trans- 
forms of this type, care must be exercised to restrict the 
range of the arguments in order that finite non-negative 
absorption is added to the problem by the transform. Other- 
wise reduced or negative "effective” total cross sections 
will result and convergence will be decelerated in the 
transform domain. 

It is possible that problems in K,\x and higher- 
dimensional geometries can be formulated in the mathemati- 
cal context displayed here for the simplest geometry. If 
transforms can be found which place an acceleration parameter 
in the proper place in the iteration matrix to reduce its 
norm, vast savings of computational effort can be expected 
for these problems. The transform must be one which leaves 
the angular discretized Sf^ equations invariant. To be of 
practical value, such a transform must not only properly 
position an acceleration parameter in the iteration matrix 
but also must not alter the source so drastically that large 
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variations occur in the <j)-domain solution with the 
accompanying unacceptable discretization error. 

Possible transform candidates are some form of the 
integrating factor for the angular discretized equations 
describing the problem after neglecting the scattering gain 
term. Transforms are not restricted to this class, however, 
and will prove useful if they meet the criterion described 
above . 

It is doubtful that a transform can be found which 
is useful for all problems in a particular geometry. 
Experience of this investigator indicates that it is more 
likely that each problem in the geometry has a "best" trans- 
form which can be prescribed based upon the source, boundary 
conditions, and the prescriber's estimate of the form of 
the solution. 

Conclusions 

1) The transform (33) applied to the half-range slab 
problem accelerates convergence of the Sj^ algorithm for 
predominantly scattering media using the improved convergence 
criterion. 

2) The transform method is successful only for prob- 

2(y-l . 

lems in which A is chosen so that A < — . 

3) Significant convergence acceleration occurs for 
values of a much smaller than a*. These smaller values of a 
produce solutions with acceptable discretization errors. 
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4) A practical choice of oi for a particular problem 
must be based upon experience and the user's estimate of 
the expected solution characteristics. 



0 



CHAPTER IV 



DISCRETIZATION ERROR IMPROVEMENT USING 
SPATIAL TRANSFORMS 



Discretization error or truncation error, as it is 
frequently called, results from discretizing an equation in 
the continuous domain and comparing the discretized equation 
to the analytic equation. For the full-range Sjij approxi- 
mation of Chapter II, this process is illustrated by 
approximating the set of ordinary differential equations 



hi 



^ 1 ^ ^ ^ V 

M; ^ + O rp .(x) j- + q;(x) 

J J yi= I J 






y = 1,2,''', hi, which are valid in the spatial continuum 
0 ^ X £ L, by the spatially discretized set 



4 -) ' » \—-i ) ■ - ! ) ' h-ui.i 



( 2 ) 



which uses a centered difference approximation for the 
derivative and a simple average for ib at each mid- 

point ■ J 

Expanding Taylor series about 

the point t^^>^<^^ting appropriately, and forming 

the proper sums produces 

^ <’1* I (3) 
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where 



(4) 



and 



2 






where X^ < < X^^j 



(6) 



Substituting equations (3) and (5) for the appropri- 
ately indexed quantities in equation (2) gives 



( 7 ) 



Now evaluating if)y(x) at equations (1) and sub- 

tracting the result from equation (7) gives the discretiza- 
tion error 

' ^ 1 



If it is assumed that the solution ipjlX} to equation 
(1) is sufficiently smooth that its third derivative exists 
and is bounded, that is. 



ip'" (x) < M (9) 

i 



for all y,x, then 



TylX^^j/^) < MO(A^) + pO(A^) 



where 



p = max 






N 



n= J 



(10) 



( 11 ) 




I 
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2 

The discretization error is clearly of order A for this 

2 

method. That is, t approaches zero as A approaches zero. 

A standard technique for diminishing the discretiza- 
tion error is to choose A small enough, depending upon one's 
estimate of the variation of the solution, that the dis- 
cretization error is tolerably small. 

Transform Method 

The above technique of diminishing A to improve the 
discretization error has a practical limit. As A decreases, 
the ntimber of spatial mesh points increases as does the 
computational cost. Consequently a choice of A must be 
based on the tradeoff between accuracy and cost. 

Another method of decreasing the discretization 
error is suggested by (10) and (11) . If a transform can be 
found which renders equations (1) invariant in a transform 
domain and diminishes the second derivative of the solution 
without increasing M, the discretization error will be 
reduced. That is, if the transform domain solution of 
equation (1) has less curvature magnitude over the slab 
than the solution i|^y(x) of equation (1), the transform 
method will improve the discretization error. This method 
will be demonstrated in the next section for a simplified 
Sj^ problem which has an analytic solution. 

Simplified Problem 
The equation 

+ q(x) 



( 12 ) 
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with boundary condition 



|o ) = 1.0 



( 13 ) 



and source 




qix.) = e 



( 14 ) 



has a solution 




-o^k/Vq 

2 ^ 



( 15 ) 



where 

( 16 ) 

and ij^ is a positive real constant. 

An S|^-type iterative solution of the same problem is 
formulated by approximating the derivative by a central dif- 
ference and approximating equation (12) at a spatial mid- 
point as before to give 




( 17 ) 



with 



= J.o 



( 18 ) 



and 





(19) 
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Here 



e » 




(20a) 



and 



d 



+ 

TE 




(20b) 



Observe that this algorithm is precisely the same form as 
the Sjtj algorithm along a |yy| ray with the exception of the 
quadrature approximation of the scattering gain term. 



The Exponential Transform 
The transform 

- otx 

4) (x) = (}) ( x) e , 



(21) 



with a > 0, can be applied to equation (12) to give the 
transform domain equation 



^0 ^dx^^ (x) + q (x) . 



( 22 ) 



Discretizing spatially, as before, gives 



.u+n _ e , f,. 






aI-^) 

<\>l 






(23) 



where 



e = 



2\x^-Aa +Aay^ 

u 



(24a) 



and 



d = 



2y^+Aa -Aay^ 

Tl 



(24b) 
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Observe that this transform places the parsimeter oi into 

Q « 1 

such a position that ^ is increased and ^ is increased. As 
shown in Chapter II, this increases the norm of the iteration 
matrix and decelerates convergence of the algorithm. Con- 
sequently the method is restricted to those problems for 
predominantly absorbing media in which the iterative 
efficiency is not so sensitive to the norm of the iteration 
matrix. For these problems, the conventional convergence 
criterion 



(X+l) 

k 

U1 




< e 



(25) 



can be used with negligible loss in iterative precision. 

The exact <j>-domain solution of equation (22) with 

source 



/ 1 ~a^x+ax x{a-o^) 

q[x) Z = e ° ^ = e 



(26) 



and boundary condition 



<j)(o) = 1.0 



(27) 



IS 



0 (x) = 



-X (Cf a) a^-y a*- / 

e ' 0 ^ 

+ g 



(o^^-ay^ ) 



(a^-y^a^^) lo^-y^a^) 



( 28 ) 



Scirnple Problem 

A problem with the parameters listed below was 
treated by the methods described above. 
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Slab width = 10.0 cm. 

R = 30 spatial intervals. 

- 1.0 cm“^. 

= 0.2 cm“^. 

= 0.9. 

e = 10"^. 

The exact solution of equation (12) for these 
parameters is shown in Figure 4.1 and is 

iplx) = . (29) 

The exact <j)-domain solutions for various values of a are 
shown in Figures 4.2 through 4.6. For values of a = 0.5 
through 0.9, the exact <J)-domain solution has much less 
curvature magnitude over the slab width than does the exact 
i|)-domain solution. Consequently it is expected that the 
transform domain iterative algorithm (23) solution will 
possess smaller local discretization errors than the 
unmodified algorithm (17) . 

Algorithms (17) and (23) were run on an IBM 360/67 
computer to verify the expected discretization error 
improvement using the transform method. Table VIII sxim- 
marizes the results. 

The fractional errors were calculated from 

Inexact _ ^U) I (30) 



exact 
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FIGURE 4.1 

EXACT SOLUTION FOR UNMODIFIED 
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FIGURE 4.3 

EXACT SOLUTION FOR a = 0.7. 




FIGURE 4.4 

EXACT SOLUTION FOR a = 0.9 
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FIGURE 4.5 

EXACT SOLUTION FOR a 



1. 1 




X,CM. 



FIGURE 4.6 

EXACT SOLUTION FOR a = 1.3 



8 




X,CM. 



DISCRETIZATION ERROR IMPROVEMENT 
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for the unmodified algorithm and 



u 



exact 

k 




<!> 



exact 

k 



( 31 ) 



for the transform domain algorithm. The errors compared, 
however, are for values of the unmodified Sj^ if^-domain solu- 
tion and the transform method i|;-domain solution. These 
fractional errors are attributed primarily to discretization 
error because the unit roundoff error for the double pre- 
cision mode of the computer is 10“^^ and the iterative 
error should be less than lO”^ due to the convergence 
criterion (25) used. Certainly fractional errors greater 
than 10”^, as most unmodified Sf^j errors are, can be attri- 
buted to discretization error. 



Discussion of Results 

The transform method produced steady discretization 
error improvement with a for all space points up to a = 0.9. 
This error improvement is due to the general decrease in 
curvature magnitude of the exact (|)-domain solution compared 
to the exact if^-domain solution over the slab width. The 
transform method can be used to diminish p in equation (11) , 
thereby reducing the discretization error. 

The value = i.i produced error improvement over the 
unmodified Snj solution but less improvement than did a = 0.9. 
Values of a > l.l produced larger discretization errors than 
did unmodified in some regions of the slab. Observe that 



r 



I 




- 100 - 



the solution <|>(x). Figure 4.4, has the least curvature 
magnitude at u = 0.9. Also at a > — 0.9, equation 

(28) reveals that (|)(x) has a growing exponential in the 
second term which dominates the solution. This suggests 
that predicting an optimum value for a requires knowledge 
of the exact solution. Obviously this is not possible in 
more complicated practical problems but some guidance is 
available from the simplified problem. 



Strategy for Determining Optimum a 

A practical value for a may be obtained from the fol- 
lowing argument. The method is only useful for predominantly 
absorbing media. These problems are characterized by solu- 
tions which "follow the source." That is, if the source 
distribution is a decaying exponential, the exact solution 
will have the form of a decaying exponential. Consequently 
a transform which changes the source to nearly a constant 
in the (j)-domain would probably produce a more slowly varying 
4>-domain solution with accompanying smaller discretization 
error. For the problem treated here, that would be 



(? ( x) e 



ax 




coyut, 



when 

a = = 1.0 . 

This is very close to the value which was found nearly 

optimum by knowledge of the exact solution. In the absence 
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of knowledge of <}> ( x ) , the above strategy based on source 
shape could be useful. Of course if partial insight into 
the expected solution is available from analytic solutions 
to similar problems or if the user's intuition is sharp, a 
value of a based on such considerations could prove fruit- 
ful. 

Application to 

The Sfj algorithm for the full-range slab problem, 
derived in Chapter II, is repeated here 




(32a) 



for Uj > 0 and 




(32b) 



for Uj < 0 
where 







(32c) 



TK 



and 





(32d) 



2A 



The transform 



= 4>y(x)e"^^ 



( 21 ) 



with a > 0 applied to equation (1) gives 
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d<P;lx) . ^ ax 

,j. -4^- + (o^-aii .(x) = f- I + o.(x)e (33a) 

J J J ^ yi= 1 J 



for iiy > 0 and 



d(J):(x) ^4 W 

~ I '^y I • -fc ~ ^ Z ^„(x)+q .(X)e (33^) 

J ■' J ‘ / J 



for viy < 0. Observe that the transform adds absorption to 
the My < 0 equations but diminishes absorption in the 
My > 0 equations. 

Spatial discretization is conducted in the manner 
previously described and coefficients of the fluxes are col- 
lected to give the Sf^ algorithm, 



\*t,j ■■ *T:\r 

^ ^ [ n= 1 



1 / 2 \ 



(34a) 



for My > 0 and 



,U*U iul-i*') J jo* T J-il *„ /’^k*UA 

*k,j ' Tj\l- ^y*k*l/2,n ^k*U2,j 

J J [ n= I J 



(34b) 



for My < 0 where 



ey - < 



a I u • I 

+ 7T^ 



a I M ; I 



for ]i J > 0 



for ^ 



(35a) 



(35b) 



and 



a I U ; I 



d; = < 



dy - — ^ for My > 0 



a I M . I 

dy + for My < 0 



(36a) 



(36b) 
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e ; ; 

Observe that a increases -r^ and t— for y • > 0. 

Tj 37 s 

Consequently the part of {V-E) 'S which corresponds to the 

)] ■ > 0 equations has a positioned such that each non-zero 

^ e; 

matrix element is increased. Conversely, a diminishes ^ 

and -I— for u • < 0. This results in a decrease in the non- 

dj j 

zero matrix elements of the iteration matrix corresponding 
to the \ij < 0 equations. These properties destroy the 
simple block collapsing technique, demonstrated in Chapter 
II, used to calculate the norm of the transform domain 
iteration matrix. Some portions of the maximum row sum are 
diminished and some are increased by a. Fortunately it will 
not be necessary to calculate | | (P-E)"^S| | for the problem 
treated here since only predominantly absorbing media are 
considered. Consequently the norm will not change signifi- 
cantly with ot nor will the number of iterations to con- 
vergence . 

Sample Problem 

Equations (34) , subject to the boundary conditions 



<j) , . = 1.0 for y . > 0 

1.1 i 



(37a) 



and 



()) . = 0 for y • < 0 



(37b) 



and source 



<?y(x) = 0 for all y. 



(38) 
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were solved on an IBM 360/67 computer. The pointwise con- 
vergence criterion 



U^l) 



- <t> 






TTF 



< e 



(39) 



was used. The problem parameters are: 

Slab width = 10.0 cm. 

R = 30 spatial intervals 
= 1.0 cm“l. 

= 0.2 cm”^. 

Quadrature set = six point Gauss-Legendre 
a = Various values from 0.1 to 1.0. 

The resulting i|;-domain solution will be compared to two 
unmodified Sjy solutions of the same problem, described by 
equations (32), but with different spatial mesh sizes. One 
used R = 600, the other R = 30. Appendix E contains the 
programs and sample output. 



Unmodified S|^ Results 

Comparison of the unmodified Si^ results for R = 600 
and R = 30 shows that all R = 600 values are either equal to 
or greater than (at least to five significant figures) the 
R = 30 solution values. Since the R = 600 solution has less 
discretization error than the R = 30 solution, the above 
fact simplifies the analysis of the effect of the transform 
method on improving discretization error. If the transform 
method i|^-domain solution values are between the unmodified 
R = 30 and R = 600 values, the method has improved the 
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discretization error. The above criterion is used to judge 
the effectiveness of the transform method in improving 
discretization error. 

Transform Method Results 

The data show that the transform method decreases 
the discretization error for the problem. Discretization 
error improvement is not pointwise uniform but, in general, 
transform method solution values approach the R = 600 solu- 
tion values as a is increased. The bulk solution differ- 
ences, displayed and defined in Table IX, diminish for each 

p ; component of the solution as a increases. The total 
J 

bulk difference also diminishes as a increases. This means 
that some integral measure of the discretization error 
diminishes as a increases through a = 1.0. A closer look 
at the pointwise data, however, reveals that at a = 1.0 
many transform method i|»-domain solution values have overshot 
the unmodified Sjy R = 600 solution values. 

Table X displays a more detailed analysis of the 
pointwise discretization error improvement achieved by the 
transform method. It shows that discretization error is 
decreased at each point for values of a up to 0.5. For 
a > 0.5, discretization error is generally improved but 
some values overshoot the R = 600 values. As a grows larger, 
more values overshoot the R = 600 values and the analysis 
loses its basis for comparison. All methods required the 
same number of iterations to converge; 10. Hence 
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of 186 have closer One value overshot R = 600 value. Agree: 

5th for y < 0 values to R = 600 generally to 2nd significant figure for 

solution than did y >0 values. Agreement generally to 1st 

a = 0.3 solution. significant figure for y < 0 values. 
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discretization error improvement was achieved at the expense 
of only R*W extra computations. 

Conclusion 

The transform method can be used to decrease the 
discretization error in Sjij problems. This method is effec- 
tive when the transform domain solution has less curvature 
magnitude over the slab width than the unmodified solu- 
tion. 

The transform used for this purpose must be tailored 
to each individual problem. The choice of a successful 
transform for a specific problem is based upon intuition 
developed for that problem by experience, knowledge of 
analytic solutions to similar problems, or guess. For 
problems in which solutions tend to "follow the source," a 
transform may be devised which converts the source to a near 
constant in the transform domain. A possible transform 
candidate is some variation of the integrating factor for 
the transport equation neglecting the scattering gain term. 



CHAPTER V 



ROUNDOFF ERROR ANALYSIS FOR ITERATIVE METHODS 
Background 

Previous chapters have discussed the concepts of 
iteration and discretization errors in the method. Each 
effect was treated with regard to its role in the con- 
vergence of the algorithm to some exact solution. The 
approximate solution was shown to approach the exact 

solution ipylx) of the angular segmented equations as 
approaches “ and A approaches zero. This concept of con- 
vergence ignores the errors that occur in any machine 
computation due to the rounding procedure employed in 
floating point arithmetic. In order to be effective, how- 
ever, an algorithm must remain immune to the accumulation of 
roundoff errors. This immunity is called numerical sta- 
bility. An otherwise convergent algorithm may be driven 
divergent if it does not possess this property. In order 
to analyze the property of numerical stability, precise 
definitions must be made. 

Well-Posed Computation 

A numerical method, which produces a solution 
approximating the exact solution of the problem, is said to 
be a well-posed computation if the solution is insensitive 
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to small changes in data or to roundoff errors. Here data 
refers to the elements of matrices involved and the Initial 
or boundary conditions and sources. Isaacson defines a 
computing problem as an algorithm or "a set of rules 
specifying the order and kind of arithmetic operations 
(i.e. rounding rules) to be used on specified data.” He 
also defines a computing problem to be well-posed if and 
only if 

i) a solution exists 

ii) the solution is unique, that is, when performed 
several times, with the same data, identical results are 
obtained 

iii) the solution depends Lipschitz continuously on 
the data with a constant that is not too large. The last 
condition requires that "small" changes in the data must 
result in only "small" changes in the solution. A well- 
posed computation possesses the property of numerical 
stability. A necessary condition that a finite difference 
scheme be convergent is that it possess numerical stability. 

Condition Number 

Roundoff errors will in general be introduced in the 
course of carrying out the arithmetic required to solve 
the system of equations 



Ax = £ . (1) 

But if the algorithm is well-posed, these errors can be kept 
within reasonable bounds. The matrix A is said to be 



If lift ntfwiid'llBiil^ 
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"well-conditioned" or "ill-conditioned" if the computation 
is or is not, respectively, well-posed. The condition 
number, v{A), for A serves as a measure of ill-conditioning. 
The condition number is defined by 

v(A) = ||A1|-||A-J|| . (2) 

Some important properties of the condition number are: 

i) v(A) >. 1 with equality only if A is orthogonal, 

ii) V (A) = V (A" M • 

iii) V (“A ) = V ( A ) . 

Suppose that the data A and £_ of equation (1) have 
small perturbations or uncertainties 6A and 6 _^ respectively 
due to roundoff error. The computer is solving the slightly 
perturbed system 



(A + 6A) (X + ^) = ^ 

instead of system (1) . Isaacson^^ shows that if 



16A|| < I/IIA-Ml , 



(3) 



(4) 



the resulting relative uncertainty in the solution, X, is 



lUM - 





MAN J ^ 



(5) 



Observe that if v(A) is relatively small, A is well- 
conditioned and "small" relative uncertainties in the data 
produce "small" relative uncertainties in the solution. 
Conversely if v(A) is very large, A is ill-conditioned and 
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"small" relative uncertainties in the data can destroy the 
solution . 

Forsythe^'* points out a popular misconception that 
the smallness of de-i(A) causes the ill-condition of A. 

That is, the condition number is not a measure of how nearly 
A is singular, hence difficult to invert. He stresses that 
the size of v( A) is a far more important criterion of the 
"badness" of computational system (1) than either the 
smallness of det(A) or the largeness of the order n of the 
system. 

53 

A theorem from Isaacson, which will prove useful 
later, is stated here. 



Theorem 2 

If the matrix A has | |A| | < 1 , then (I ± A) is non- 
singular and 



;HIa 



1 M (I t A]'1\\ < 



- 1 



All 



(6) 



Application to Iterative Methods 

Isaacson, referring to iterative methods for 
solving systems of equations (1), states "One of the intrin- 
sic advantages of such methods is the fact that errors, due 
to roundoff or even blunders, may be damped out as the 
iterative procedure continues". The following analysis 
provides an insight into the mechanism of this intrinsic 
advantage . 
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The interative S/^ algorithm, described in Chapter IT, 
has a matrix representation 



^ . (7) 

At each iterative step, the data D, E, S, and are 

known and the machine is expected to invert (D-E) exactly 
and multiply the result onto the right-hand-side of (7) , 
that is, solve for 



= (D-E) ■ + (D-E)'^£ 



( 8 ) 



precisely. The computer, however, has a limited storage 
capacity for each digit and must round off the input data 
to the last significant digit of its capacity. Consequently 
the uncertainty in the stationary data is 6D, 6E, 6S, and ^ 
which are of the order of the unit roundoff error of the 
machine (10“^^ for the IBM 360/67 in double precision mode). 
These uncertainties are propagated through the arithmetic 
process of inverting (D-E) and matrix multiplying the 
result onto appropriate right hand side vectors. Conse- 
quently even the input data has its uncertainty. Sip , 

carried with it in equation (7) . Therefore the exact system 
to be solved at the iterate is not (7) but 



(D-E)]jj^ 






(9) 



where the previously computed solution, r is 






( 10 ) 
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But now the machine must solve the perturbed system 

(D+6D-E-6E) = (S+6S) + (11) 

where is the roundoff error resulting from machine 



handling of the data ^ 



l-i) 



It is the same order of magni- 



tude as 6 d, 6e, 6S, and These uncertainties in the input 

data produce an uncertainty ^ ^ in the solution. Using 

equation (11) , it is possible to derive an upper bound on 
the norm of the relative uncertainty in the solution 



I 

TTTiMTjJ 



at each iterative step. 



First Iteration 

Suppose the initial guess for algorithm (7) is 

= 0 ( 12 ) 

with no machine error. That is, the first machine iteration 
is on The first iterate is 

(D-E)^^ = g_ (13) 

but the computer solves the perturbed system 

(D+6D-E-6E) (^f ^ ^ * ) = £_ + ^ , (14) 

Factoring out (D-E) from the left-hand-side and multiplying 
by (D-E)~^ gives 

[l+ (D-E) ■ ^ (6D-6E)] = (D-E)"M2^+^) . (15) 
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Recognizing that 

= (D-E)’^£ 



(16) 



from (13) , cancelling this equality in (15) and collecting 
terms in on the left-hand-side produces 

[l+ (D-E) " ^ (6d-6E)] ^ ^ = (D-E)'^ (6E-6D)i|^^^^-»-(D-E)~^6q . (17) 

Hence 

^ [l+(D-E) ■ ^ (6D-6E)] 

But due to Theorem 2, 



- 1 



j(D-E) (6E-6D)i|^^ ^ ^ + (D-E) " ^ §j^ 



(18) 



I I [l+ (D-E)"^ (6D-6E)] I I 1 



I- I I (D-E) ■ ' (6d-6e) 



(19) 



subject to the restriction that 



I 1 6D-6E| I < 



I I (D-E)-' I I 



( 20 ) 



Hence 



(-||(D“B)''{«d-«e)|| 



|l|(D-E)''|l'll«D-«E|Mli‘''|| + I|(D-E)*'||-Il«lll| 



Now divide by 




and multiply by 



1 


(D-E) 




1 


(D-E) 





to obtain 



< V I ,n, 

J- I I (D-E)"N6D-6E) I ill I (D-E) I I I I (D-E) I I • I li*’h iJ * 



where v is the condition n\imber of the matrix (D-E) , that 
is, 
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V 



||(D-E)||.||(D-E)-'|| . 



( 22 ) 



But from (13 ) , 



Hill < ii(D-E)i|.|ii‘'’ii 



or 



< 



( 23 ) 



I|D-E||-Ill'"ll Hill 



From restriction (20) and the fact that 

I I (D-E) ■ ^ (6D-6E ) I I < II (D-E) " M I * I I 6D-6E | | < ? , 
r- I I (D-E) ■ ^ (6d-6e) I I > I- I I (D-E) “ ^ I I • I I 6D-6e I I hence 



r- I I (D-E) ■ ^ (6D-6E) II I- I I (D-E) ■^11*11 6D-<SE| | 
Using (23) and (24) , inequality (21) becomes 




< 



J-v 




(25) 



Second Iteration 



The second iteration treats the problem 




(26) 



but the machine solves the perturbed system 



(D+6D-E-6E) = (S+6S) ^+611;^ ^ ) + (^+^) . 



(27) 
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Following the same procedure used during the first itera- 
tion, 

V || |6E-iP|| ^ ll^.s-11 I l ie" II , l|S>*-’i|l lll<i|l 




A bound on the quantity 
That is, 



Ill's’ll 



can be derived from (26) 



(21 - 7 ( n - 7 

= (D-E) S± + (D-E) 1 



or 






[(D- 



E)'^s] - (D-E) V] • 



Consequently 

IW ” 



< I I |(D-E)'^s] ^1 I [| 



(D-E) 



- 7 



or 



t 



(7) 






(D-E) "'sj 



- 7 



, ^ l l(D-E)'MM I i U 



(2 



(29) 



A bound on ; ^ r-; on the R.H.S. of (29) can be derived by 

lit' 'll 

using the unperturbed system 



= (D-E) ■ ^ ^ + (D-E)~^q 



(30) 



From equation (13) 






(7) 



(D-E) ■ 



(31) 
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which substituted into (30) gives 

[l + (D-E) " ^ ^ ' 



or 



= [l + (D-E)‘^5]"^^(^) . (32) 



Using Theorem 2, 

|uin|| < Hit's'll - 

subject to the restriction that 

I I (D-E) " ^5 I I < / . (33) 



But recall that (33) is the condition that the Sj^ algorithm 
must meet to guarantee convergence so it is met a priori. 
Consequently 






But from (23) , 




I 


D-E 


1 




kl 





so the above inequality becomes 



? ^ I I D-E I I ^ (34) 

lli''’ll ' llll|-f/-||(D-E)-'s||) 



Using (34) on the R.H.S. of (29) produces 




;+v- 1 


(D-E) 




r- 1 1 (D-E) si 1 





( 35 ) 
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Also using (34) for the factor 
the bound 



Jl, 

'P 



n 



2 } 



in (28) provides 



nay’ll ^ iiD-Eii-iia^”ii 

Ill's’ll " ||4||-(I-||(D-E)-'S||) 



Now applying (34 ) , (35) , and (36) to (28) produces 



I'lpl 1 1 ' 





1 1 6E-(SD| 


[1 , ' 


IIS5II .p(l»v- 1 l(D-E) ‘ 'S| n 


(-'«• I'iAd-AeII 1 


1 Id-eI I 


[l-l 1 (D-E)-'s| 1] 


||d-e|| ||(d-e)''s|| 


itb-El f 









il 



( 37 ) 



where 



p = I I (D-E)'^S| I . I I [(D-E)'^S] ’^1 I 



is the condition number of the iteration matrix. 



(38) 



The General Iterate 

Suppose the iteration has been completed. That 

has been computed. A bound on the 
relative uncertainty in the iterate solution can be 

derived from recognizing that instead of solving the system 

(D-E)tli^'^^^ * = ^ (39) 

the machine solves the perturbed system 

(D+6D-E-6E) = + + . (40) 

Proceeding in the same fashion as in the previous section. 
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, V IIISE-6DII . lUSlI 

■ T: ^-T||^grp- 1 ii^Bji ■ twT 



I 1 1 

Ili'-'MI 



S«5S| I . 
D-E I 



IIMII 






( 41 ) 



From equation (39) , 

ij) (-<-+?) = (D-E) ■ + (D-E)"^£ 



or 



= [(D-E)-'sJ‘^[^^f^^^J - (D-E)-^£j 



Consequently 

Hi'"’ 1 1 



•c 

TTHT 



F/ 




1 1 (D-E) ” ^ 1 1 - 


‘iiiir 


(D-E) Sj 


• 


r “ 07 ^ 





(42) 



The factor 



i 



imj 



on the R.H.S. of (42) can be 



bounded by recalling that the unperturbed system obeys 



. f a’”(d-e)-'£ 

m-0 

where for convenience of notation 



A = (D-E)"^S 



(43) 



( 44 ) 



Consequently 

2. • (D-E) [l ♦ A » ♦• • •+ A-^J 



-'^u*n 
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But if condition (33) is met, 



(I-A)'^ = I A"' 
m-0 



so 

[l+A+A^+- • • +A-<^] = fd-A)"' - A^'^' - ]'^ 

= [(I-A)"' - A^’'^ (I + A + A^+* • ^ 

= [(l-A'^^') (I-A)-^]'^ 

- (T-A)(I-A^^^)-' . 

Hence 

q = (D-E) (I-A) . 

Consequently 

11,11 < IId-eII-(i^IIaII)-II»i^*’'II 

’ - lU'"' II 

due to Theorem 2 if 

IIa'""|I < ' 

which is certainly true if condition (33) is met. Since 

IIa^'MI < IIAII'^*' < I , 



(45) 



)-||A^*'|| i >-||A||^*' 
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and 



1 ? 

i-l|A-i*'|l - i-I|a|F^ • 

Using (46)/ inequality (45) becomes 

1 1^1 1 , l|D-E||-(r>||A||)-||tl^''l|l 

1-llAll'"*' 



( 46 ) 



or in terms of the iteration matrix norm 



1 


, l|D-E| 


•[/ 


+1 1 (D-E)"^S| |] 




1 Ikll- 


I- 1 


1 (D-E) ■ h\ 1^^^ 


(47) 



Applying (47) to the R.H.S. of (42) gives 



nils'll p 

lll'A'Ilii ^ ||,o.E,-'s|| 



‘J + v (l+l I (D-E)''s| l) - | I (D-E)"'s| 
I-||(D-E)-'S||^*' 



(48) 



Now applying (47) and (48) to (41) produces 



111*"'” II 


V |||ae-ad|| 


I 


11^511 


1 1 (D-E) "'si i) - 1 1 (d-e) ■ 'si r"] 


llt”"’ll ", 


V 1 |6D-6e1 I ] 1 |D-e| I 1 

-tt^n i 


f- 1 1 (D-E) “ 


1 |D-E| 1 


||(D-E)’S|| 



||.S*AS||.|>M|(D-K)-'s|||.||^j'l|| 

* TmT“ 




(49) 



This expression is valid for X = 2,3,'''. 



Analysis of Error Bound 

Observe that in expression (49) , the only terms which 

are affected by the number of iterations are | | (D-E)'^S| 

and 1 . As -c increases, the former magnifies 

I - I I (D-E) - JS I ^ 
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the bound slightly and the latter diminishes the bound. The 
latter effect is predominant. Therefore the bound becomes 
tighter as iteration proceeds and roundoff error propagation 
is diminished. This effect is not strong enough to force 
the relative uncertainty to vanish, however, as -L approaches 

Certain terms and factors in (49) can be neglected 
with respect to more dominant quantities in order to obtain 
a qualitative view of the bound. That is, 

I '* 1 1 y 1 1 I«e-6d| I ^ I |6S| I ^ I |s<as| I I ImI* * 1 1 ^ I llll (50) 

~ , V I hp-6E| I [ ||D-E|| ||d-e|| ||(D-E)''S|| ||£|| Hill [ 

I |D-En 



contains the dominant quantities in (49). 

The condition of the matrix (D-E) has a dominant role 
in the roundoff error bound. If v is very large, the rela- 
tive uncertainty in the iterate solution will be significant. 



If V is large enough that the factor 



V ■ 



6D-6E 



D-E 



is close 



to unity, the relative uncertainty can be extremely large. 
For this to occur, v would have to be of the order of 10^'^ 
for the double precision mode, extremely ill-conditioned. 
This event is unlikely. Since v, which is always greater 
than unity, appears to the second power in (50) , it is the 
dominant quantity in the bound. 

The condition number of the iteration matrix, p, also 
plays an important part in the roundoff error bound. Its 
effect is subordinate to that of v since p appears only to 
the first power. But if (D-E)”^S and (D-E) are both 
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ill-conditioned , the quantity v^p could cause significant 
roundoff error propagation in the solution. 

The error bounds derived here display another 
important effect. Expressions (25), (37), and (49) each 

contain the uncertainty ||6D-6 e|| in important positions. 
This can cause the exact cancellation of roundoff errors 
during the iterative process. That is, | |6D-6 e| | vanishes. 
This gives vivid corroboration to Isaacson's observation 
quoted previously. The presence of | |6D-6E| | in the 



dominating factor 



V 



1 



V • 



6d-6e 



reduces it to v when exact 



D-E 



cancellation occurs. This effect is also felt through the 

term I |6E-6 d| | . 

I I D-E I I 



Observe that each term in the braces of (49) has a 
multiplier which is extremely small. That is 6E, 6D, 6S, 
6q , and ^ are all of the order of the machine unit 

roundoff error. The factor | | ^^ | | is generally on the order 
of unity. From equations (23) , (24) , (25) , and (26) of 

Chapter II, the infinity norm of S is 



4 hi 

1^11 • ^ I 



= O' 



( 51 ) 



From equations (16c), (16d), (21), and (22) of Chapter II 

the infinity-norm of (D-E) is 
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||D-E| 




= max 

j 



A 



^ I I max 

A 



( 52 ) 



As a consequence of the above facts, the quantity 






/-■ 



V • 




1 6D-6E 






|D-E| 1 



must be of the order of 10^^ for the relative uncertainty to 
be significant. Assuming that | |6D-6 e| | = 0 and p ~ v, this 
means that 



or 

V r 10 ^ 

before significant roundoff error is propagated to the 
iterate solution. These extremely large condition numbers 
simply do not exist for physical problems of interest in the 
application of the method. 

Estimating the Condition Numbers 

The condition number of the iteration matrix is not 
readily estimated. Although the infinity-norm of (D-E)"^S 
can be calculated precisely, as discussed in Chapter II, 
calculation of | | [(D-E) “ ^Sj " ^ | | is not tractable. 
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Consequently an estimate of p is not available nor is an 
upper bound on it. The following discussion assumes that p 
is small enough that its influence on the roundoff error 
bound is less than that of v. 

Estimation of v is possible. Equations (27) and (28) 
in Chapter II reveal that the infinity norm of (D-E)“^ is 




or 



II(d-e)-M| 



1 



7^ 



(53b) 



where R is the n;imber of equi-spaced spatial intervals, A, 
and 



0 < 3 




(54) 



Here 



e , 



2 I u . • - Ao' 

' '^j ' m-tn 

Ta 



(55a) 



and 



a, . 



m-cn 

~fA 



+ Ao' 



(55b) 



From (53) , 



I I (D-E)‘^l I 



/ . ) 



< 



( 56 ) 
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Consequently from (52) and (53) 



V 



^ 1 I max 1 
Ad. 1 - 6 



( 57 ) 



where is the index for the minimum absolute ordinate in 
the set {jJy}. 

A tight bound on v can be obtained by neglecting 3 • 
That is, 

^ I •^y I max 

V < ^ 

Ao^ 



Since all problems must have 



V < 



^y ' max 




< /, 



(58) 



is a useful estimator for v. For the problems considered 
in Chapter II, = 1.0 cm" ^ and A = / / 3 cm so 



V < 6 . 



(59) 



Hence (D-E) is a very well-conditioned matrix for the 
problem considered. Consequently roundoff error propagation 
was insignificant in the investigations conducted. In 
practice, roundoff error propagation has not been experi- 
enced by this investigator which renders the algorithm 
unstable. This must be a consequence of the well-conditioned 
state of the matrices (D-E)'^S and (D-E) as well as to the 
roundoff error cancellation property of the iterative method. 
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Effect of Exponential Transform 

The presentation in Chapter III reveals that the 
transform 



ip-(x) = (60) 

j J 

applied to the equations, produces an iterative algo- 
rithm with the same form as that of the unmodified 
algorithm but with a modified source and a different 
iteration matrix, [V-E)~ ^ S . As a consequence, the roundoff 
error analysis, previously presented, carries over into the 
transform domain in its entirety. That is, expressions 
(25 ) , (37 ) , and (49 ) are valid for 



1 




1 




|±(-i*n 1 







if (D-E) and (D-E) ’ are replaced by iV-E) and (f-E) 
respectively and ^ is replaced by its transformed counter- 
part. All of the above quantities are defined in Chapter 
III. 

From equations (18) and (19) of Chapter III it is 
clear that 



M(P-E)|| = M (D-E) I I 



(61) 



since terms involving ot cancel. From equations (23) and 
(24) of Chapter III, 



J-Y . 



(62) 
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where 



0 £ y 




(63) 



and 



Z I M , I . -Aa'^-Aal y . I . 

' j ' m-cn ' j ' 

f5 



(64a) 



Z I y . I . +Aa'^+Aa| y . | 



- 



j ' mA^n 






j ' m^n 



(64b) 



Observe that in problems for which the transform method is 
effective, i.e. Q.j >. 0 for all j, and 





of equation (55b) and 



(65) 



Y < 3 



( 66 ) 



of equation (54) . 
(62) , 



Since 







the predominant factor in 



|||P-E|-'|I < ||(D-E)-'| 



(67) 



Consequently the condition number of IV-E) is smaller than 
the condition number of (D-E) . As previously discussed, the 
condition number of this matrix is the predominant factor 
in the stability of the iterative algorithm. Consequently 
it is expected that the transform domain iterative algorithm 
is "better" conditioned than the unmodified algorithm. 



CHAPTER VI 



CONCLUSIONS 



The Sj^ algorithm in one-dimensional plane geometry 
has been formulated within a mathematical framework which 
allows detailed insight into its convergence properties. 
Knowledge of these properties permits the derivation of an 
improved convergence criterion which, if met, guarantees 
that the fractional iterative error is arbitrarily small. 
Utilization of this criterion. 



m 



- I I (D-E)'^s| i] 



( 1 ) 



is especially important in problems for predominantly scat- 
tering media for which | | (D-E)”^S| | is close to unity. 

The matrix formulation of the iterative algorithm 
allows calculation of the infinity-norm of the iteration 
matrix, | | (D-E)~^S| |, for homogeneous media problems in 
which scattering is isotropic in the Laboratory Coordinate 
System. Imposing the condition 

I I (D-E) ■ ^S| I < 7 (2) 

leads to sufficient conditions for which convergence of the 
S|^ algorithm is guaranteed. These conditions impose 
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constraints on the values used for the numerical quantities 

of the problem. A, ly.l , , , and o^. The constraint 

' j ' m^n 

A < * (3) 

a 



is sufficient to guarantee convergence and allows a limited 
amount of negativity in the matrix elements of the iteration 
matrix. The constraint 



A < 




m-cn 



o 



(4) 



renders the iteration matrix non-negative and is sufficient 
to guarantee convergence for all physically realizable 
problems. The solution for such a system is non-negative 
for non-negative sources. 

A spatial transform of the type 

Ip .(x) = <p .Ixj ^^(x) (5) 

J J 

was used on a special half-range problem. This transform 
places the acceleration parameter, «, in positions within 
the iteration matrix such that each non-zero matrix element 
is reduced. This diminishes the norm of the iteration 
matrix and accelerates convergence of the Sj^ algorithm. 
Reductions of almost 25% in the number of iterations to con- 
vergence have been realized for sample problems using the 
exponential transform 

ax 



ipj(x) = (j)y(x)e 



( 6 ) 
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where a > 0. A practical maximum value of ot for this 
transform was found to be 



“ ■ ■■ A|p,-1 . 

J m-cn 



(7) 



A practical useful value of « is much smaller than u* and is 
chosen as a compromise between the user's estimate of 
tolerable discretization error and the reduction in calcu- 
lational effort. It is possible that Sj^ problems in higher 
dimensional geometries can be formulated in the mathemati- 
cal context displayed in this investigation for plane 
geometry. If transforms can be found which leave the 
angular discretized equations invariant and reduce the 
norm of the iteration matrix, vast savings of computational 
effort can be realized for these problems. Possible trans- 
form candidates are some form of the integrating factor for 
the angular discretized equations describing the problem 
after neglecting the scattering gain term. 

Nothing done here suggests that a transform can be 
found which is useful for all problems in a particular 
geometry. Experience of this investigator indicates that 
it is more likely that each problem in the geometry has a 
"best" transform which can be prescribed based upon the 
source, boundary conditions, and the prescriber's estimate 
of the form of the solution. 

The spatial discretization error in the S^j algorithm 

is proportional to the second derivative of the solution 

ip.(xj. The spatial transform 
i 
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ij; . (x) = (|) . (x) e" (8) 

J 

with a > 0 was used on a one-dimensional problem to decrease 
the discretization error. Discretization errors are shown 
to decrease because the ({)-domain solution has less curvature 
magnitude than the unmodified Sjy t|;-domain solution. 

The transform used to decrease discretization errors 
must be tailored to each individual problem. The choice of 
a successful transform for a specific problem must be based 
on intuition developed for that problem by experience, 
knowledge of analytic solutions to similar problems, or 
guess. A possible candidate is some variation of the 
integrating factor for the problem. 

An upper bound was derived for the relative uncer- 
tainty in the iterate solution of the S|^ algorithm due to 
roundoff errors. This bound is expressed in terms of the 
condition numbers of the matrices (D-E) and (D-E)'^S. The 
bound reveals that roundoff error cancellation can occur in 
the iterative process due to the factor | | 5D-6E | | . If 
cancellation does not occur, however, roundoff errors can 
be magnified during the iterative process if the matrices 
(D-E)"^S and (D-E) are extremely ill-conditioned. Such an 
event is improbable. The bound on the fractional uncertainty 
for the iterate solution diminishes as iteration proceeds. 

It does not, however, vanish in the limit of an infinite 



number of iterations. 
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For the problems treated here, (D-E) is well- 
conditioned and its condition number is calculable from the 
problem parameters. The condition number of the iteration 
matrix is not calculable because | | [{D-El'^sJ ^| | is not 
calculable. In the problems treated here, propagation of 
roundoff error was not significant. 
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APPENDIX B 



Properties of ||(P-E)'^S|| 

The norm of the transform domain iteration matrix is 
described by 




in which 



0 <_ Y = 




1 



and 



2 ( VI • ) - -Aa'^-Aa ( y • ) - 

/C (Ti/cn X. 

2A 



- 



2 ( y . ) , +Aa^ + Aa (y - ) - 

^ m^n 'C m^n 

?A 



( 2 ) 



(3) 



(4) 



and R is the number of spatial intervals in the mesh. 

The fact that | | ( t>- E ) " | | is a monotone decreasing 

function of ot is verified by analysis of three regions of 
a; small oi near zero, at « = «* , and for a > a* where 



a 



* _ 



2 (y .) . -Aa^ 

A. mA,n 

2 (y ' I - 

mA.n 



(5) 
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The derivative of | | ( t>- E ) ~ | | with respect to a is 



d||(P-E)-^S|[ 
da 



A. 



(1-Y^\ R-I 

(T^ )-y 



4 2 






(J-Y) 



+< 



- (R- ?)y 



4' 



for , fg V 

2y .-Aa^-Aay .>0 . 






/ (6b) 

2y ■ -Aa^-Aay • = 0 . 



The relationships 



d(rf^) _ ^ 

da " 2 



(7) 



^ = 
da 



Ad . 

■t 



(8) 



for 2y.-Aa -Aay. > 0 have been used in the above determi- 

A, A, 



nations . 



When a is small and 2y--Aa -Aay; > 0, 



l^-tl 2u^-Ao^-AaU^ 
y = = = 1-e 

2y; + Aa +Aay . 

A^ A, 



(9) 



where 



e = 






) + 



2 I 



( 10 ) 
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Expanding for small e. 



“ 1-Re 

and substituting into (6a) produces 

i 

o u 



( 11 ) 



da 






* 


A 2 
a y . 


- 


|2R-Ike|R-n 


- — T (R-n 
2Ul 


R|2*e)-(f*e) 



(12) 



This quantity is negative for all R. When R is large, as 
in most practical problems, | | (P-E)”^S| | has a large nega- 
tive slope. 

t 

When a = a*, 2vi;-Aa -Aay- = 0 and Z: = y = 0. In 
this case 



d| I (P-E)'^Sl 1 

da ■ ■ ~T~ 



(13) 



from (6b). Observe that the slope of | | (P-E)“^S| | is still 
negative but is much smaller in magnitude than when oi is 
small near zero. 

By definition (5), a* is the value of a when Z: = 0, 
that is. 



Aa* (y>) ' +Aa'^-2|y-) - = 0 

Wycn m^n 



(14) 



So at a*, Y = and 

For a > a*, let 

a = a* + e 



■4 A •* 

a _ Aa 

31 ■ 2 



(15) 



(16) 
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where 0 < e < “. Then 



Y 



Ae 

4 + Ae 



(17) 



I 2A / / \ 

" nri — (18) 



and 



(P-E)’^S| I 



Ag^ 

^ ^ ^ m-tnl 



(Ae) ->• 2 (Ae) ^ 
( 4+Ae ) ^ 



(19) 



The derivative of | | (P-E)"^S| | with respect to Ae is, 



d| I [ V-£)~h 
dTAii — 



2(2R-J) (Ae)*^"^ + S { R - 1 ] 

RTl 

(4+Ae ) 



( 20 ) 



which is negative for all e, 0 < e < “. 

Equations (12) , (13) , and (20) reveal that the norm 

of the transform domain iteration has a negative slope for 
0 < a < “>, consequently it is a monotone decreasing function 
of a. 

Equation (19) displays that 



Um \ \ ( V - E )' h \\ 
a -+«> 



Ae-*-°o 



(ac" 


(Ae)'^ + 2 ( Ae )^'^ 




12(m-) . 

[_ n \ A.n 


' V 

(4+Ae) 


j 



since 
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t-im 

Ae->-“ 



(Ae)^ 2(Ae) '^~^ 

(4 + Ae)'^ 



R/ 

ITT 



by R applications of -C' Hospital's rule. 
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